Click here to Skip to main content
Click here to Skip to main content

Elliptic integrals

, 5 Apr 2013 CPOL
Rate this:
Please Sign up or sign in to vote.
Description and code for calculating elliptic integrals



Introduction 

Elliptic integrals of the first and second type is given in code, with some example of their implementation. I had recently a need for these types of integrals, both complete and incomplete, and such implementations was not easy to find (I couldn't find anything in C# or VB, as usual, and the C++ versions I found didn't work properly).

As I was searching the web for solutions I realized that there were few articles that dealt with the implementation of the incomplete integrals and how you could use them to solve some specific problems. So I also included the basic theory behind the implementation.

Background   

The definition of an ellips stems from the Greek word elleipsis that has the meaning “an oval”. It reffers to the shape of an object formed as shown below.

The interesting thing about an ellipsis is that in many ways it is  formed by a deformed circle, and on further investigation some formulas are actually just an expansion of the circular formula. The different formulations for forming an ellipsis and a circle is in Cartesian coordinates quite similar:

ShapeCartesian generation function
Circle(x2/R2) + (y2/R2) = 1 
Ellipse(x2/A2) + (y2/B2) = 1 

In the table R is the radius of the circle, while A and B are the two different "radius" of the ellipse. If one looks at the two equations one might expect some similarities in calculating its properties. We start off by the area calculations:

Shape  Area 
Circle Pi*R2 
Ellipse Pi*A*B  

This seems like a very natural generalization considering the formula of the circle. We proceed with volume calculations, and this time the names of the circle and the ellipses shifts to an Ellipsoid and a Sphere:

Shape Volume
Sphere (4/3)*Pi*R3
Ellipsoid (4/3)*Pi*A*B*C

Now, there is one area where the formulas of circles and spheres, does not expand to Ellipsis and Ellipsoids, and the is concerning the arc length and surface area respectively. (I have left out the example of the spheroid, were the surface area could be found using trigonometric functions.)

The first on to pose a solution for this problem was Wallis. However due to the nature of the progression I will first cover another basic tool in mathematics, namely integration. 

Integration

This might surprise people that have learned how to integrate functions in school, but the inventor of the formula for integration is actually attributed to Fermat (and a different version by Wallis) around 1640, thirty years before Newton and Leibniz discovered differential calculus.

What Fermat discovered was a way of calculating the area below a hyperbolic curve. He simply derived the cure by using a power series. What he did was this:  

The picture is taken from the article by Joel Rosenfeld and can be viewed here (a slightly different approach of showing the same is given in this article. The latter follows Archimedes approach of exhaustion). The function is y = x k , and with choosing r < 1 one can form a geometric series were the choice of r can lead to infinitesimal  summations given that r approaches, but are not equal to 1. Each rectangle have an area given by the formula below:

And if you follow the further steps in the aforementioned article, you would eventually get the formula for the integral from 0 to a (meaning the area under the curve): 

Both Fermat and Wallis missed the fact that integration and differentials is the except opposite of each other, and that was indeed the discovery of Newton and Leibniz.  For a more complete history behind the integration you could see the article here. It was precisely Wallis's work that lead him to be able to calculate an approximate formula of the arc length of the ellipse.

Fermat's proof is actually vital for several problems in engineering, as one is now able to solve a multitude of problems involving power series and expansions of series using Taylor or MacLaurin series expansions on problems. It is actually this very technique that is used to find the arc length on elliptic shapes. 

Complete solution to the pendulum equation - Elliptic Integrals of the First kind 

This type of equation was not the first kind to be derived, despite its name, but it is possible to write the exact solution to the pendulum equation using Elliptic integrals of the First kind. The derivation is quite long, and I'm just going to show you the main steps to follow to understand how this is. The complete derivation of the equation could be viewed in the video series that is given in the references.

 

If we have a pendulum that is held by a string and with a maximum angle of displacement, we can write the equation as follows: 


After this first section, most textbooks now exchange the sin x with x as an approximation. This works well for small angles as the sin function could be expressed in a series x were x is the first coefficient.  However if one does not want to do that then one can reformulate the problem by integrating the last previous equations, and use the start conditions to find the integration constant:

 

By using the definitions we have now arrived at an elliptic kind of integral. But to get it on a standardized form on can use the following trigonometric substitution: 

And the resulting equation becomes:

 

By changing the integrand one can transform the above equation (incomplete elliptic function of the first kind, we call it incomplete when it is not 90 degrees) to the normalized complete elliptic integral of the first kind, with the constant k defined as:

And the resulting expression with the complete elliptic integral of the first kind on the left side is:

 

The derivation of the expression is quite sketchy, and I have deliberatly skipped some steps. The point here is simply show one of the uses of the elliptic integral of the first kind. The Complete formulation of the time perios T can therefore be calculated using the expression above, and the complete integral of the first kind by the letter E (the incomplete ellptic integral of the first kind is usually detonated by the letter F):

This result could now be compared with the approximation result that is usually done in most textbooks with a small angle, meaning much less than radian angle 1:

The above equation was first found by Christian Huygens, who also proposed a more complex pendulum in clocks. The reason is that a pendulum would, due to friction and other forces in the real world, slightly change the period time, and you don't want that in something as a clock. He found that in order to have the time period independent of the angle, you would have to move the pendulum in an cycloid orbit, instead of a circular one. Other people later found way of achieving this kind of orbit in a practical way.  

The actual code for the elliptic integral can be found by the use of the MacLaurin series expansion, and the code for the complete elliptic integral of the first kind is: 

    ''' <summary>
    ''' Return the Complete Elliptic integral of the 1st kind
    ''' </summary>
    ''' <param name="k">K(k^2) absolute value has to be below 1</param>
    ''' <returns></returns>
    ''' <remarks>Abramowitz and Stegun p.591, formula 17.3.11</remarks>
    Public Function CompleteEllipticIntegralFirstKind(ByVal k As Double) As Double
        Dim sum, term, above, below As Double
        sum = 1
        term = 1

        above = 1
        below = 2

        For i As Integer = 1 To 100
            term *= above / below
            sum += Math.Pow(k, i) * Math.Pow(term, 2)
            above += 2
            below += 2
        Next
        sum *= 0.5 * Math.PI
        Return sum
    End Function

This is standard textbook solutions, however most solutions are restricted to only small displacements of theta, since it is very easy to solve this equation. If you instead start to develop the solution for all angles up to 90 degrees (meaning that the pendulum string is parallel to the ground) you will have to resort to the Elliptic integrals of first kind. If you don't you would have some rather large errors in your solution, as can be seen here as a function of angle.) 

The equation of the pendulum is actually quite vital in science like acoustics, were the speed of sound could be derived from a pendulum equation (This was actually Newton that first postulated, but it was heavily disputed by Lagrange, who found a more complicated way of calculating the speed of sound. To Lagrange's amazement the solutions were however exactly the same. For the full story you could read the introduction in Lord Rayleigh classics Theory of sound Volume 1).

The elliptic arc length - Elliptic Integrals of the second kind

Elliptic functions first appeared in 1655 when John Wallis tried to find the arc length of an ellipse, however elliptic integrals got its name from Legrendre based on the fact that Elliptic integrals of the second type yields the arc length of an ellipse.   

We start off by a definition of an ellipse which in Cartesian coordinates yields (that would give you any corresponding x or y value given one input given that a > b):   

      

The equation above is in the Cartesian coordinate system (x,y,z), and to calculate the arc length we would have to transform the problem into a Polar coordinate system (r,theta). This was indeed the system used by Gauss, Lagrange , Abel and Jacobi, as they all discovered that you could write the function of an ellipsis as two separate circles, one for each length from origin of the ellipse. What I mean is this:

Here the small r are the radius of the small (red) circle, that represents the shortest width of the ellipse, and the  capital letter R represents the long width of the ellipse in blue. Each point on the ellipse could then be found using the different side lengths.   

The reasons for doing this is that we now can form an equation that describes the change of the ellipse based on the angle. However we have two different circles, which we need to add up using Pythagorean theorem. We essentially get this equation after some manipulating :

 

The arc length change per angle is now found, it just remains to integrat over all possible angles. However, the complete elliptic integral is taken from 0 to 90 degrees, but since that is exactly 1/4 over the total arc length due to symmetry. The resulting equation is therefore:

But it was not until the late 1700’s that Legendre began to use elliptic functions for problems such as the movement of a simple pendulum and the deflection of a thin elastic bar that these types of problems could be defined in terms of simple functions. 

    ''' <summary>
    ''' Return the Complete Elliptic integral of the 2nd kind
    ''' </summary>
    ''' <param name="k">E(k^2) absolute value has to be below 1</param>
    ''' <returns></returns>
    ''' <remarks>Abramowitz and Stegun p.591, formula 17.3.12</remarks>
    Public Function CompleteEllipticIntegralSecondKind(ByVal k As Double) As Double
        Dim sum, term, above, below As Double
        sum = 1
        term = 1

        above = 1
        below = 2

        For i As Integer = 1 To 100
            term *= above / below
            sum -= Math.Pow(k, i) * Math.Pow(term, 2) / above
            above += 2
            below += 2
        Next
        sum *= 0.5 * Math.PI
        Return sum
    End Function
 

The imcomplete Elliptic integrals of the first and second kind 

The approximation which are shown for the complete elliptic integrals above is usually good enough for most purposes, but when it comes to the incomplete versions it is advised to use a more modern approach. That is namely the Carlson symmetric form, and to do the first and second kind one only needs to implement Rf and Rd of those two. 

The main article were most of the information regarding the implementation comes from the new edition of "Handbook of mathematical function", which is available online here.

The implementation of Rf is given below: 

    ''' <summary>
    ''' Computes the R_F from Carlson symmetric form
    ''' </summary>
    ''' <param name="X"></param>
    ''' <param name="y"></param>
    ''' <param name="z"></param>
    ''' <returns></returns>
    ''' <remarks>http://en.wikipedia.org/wiki/Carlson_symmetric_form#Series_Expansion</remarks>
    Private Function RF(ByVal X As Double, ByVal Y As Double, ByVal Z As Double) As Double
        Dim result As Double
        Dim A, lamda, dx, dy, dz As Double
        Dim MinError As Double = 0.0000001

        Do
            lamda = Math.Sqrt(X * Y) + Math.Sqrt(Y * Z) + Math.Sqrt(Z * X)

            X = (X + lamda) / 4
            Y = (Y + lamda) / 4
            Z = (Z + lamda) / 4

            A = (X + Y + Z) / 3

            dx = (1 - X / A)
            dy = (1 - Y / A)
            dz = (1 - Z / A)

        Loop While Math.Max(Math.Max(Math.Abs(dx), Math.Abs(dy)), Math.Abs(dz)) > MinError

        Dim E2, E3 As Double
        E2 = dx * dy + dy * dz + dz * dx
        E3 = dy * dx * dz

        'http://dlmf.nist.gov/19.36#E1
        result = 1 - (1 / 10) * E2 + (1 / 14) * E3 + (1 / 24) * E2 ^ 2 - (3 / 44) * E2 * E3 - (5 / 208) * E2 ^ 3 + (3 / 104) * E3 ^ 2 + (1 / 16) * E2 ^ 2 * E3

        result *= (1 / Math.Sqrt(A))
        Return result

    End Function

And for Rd: 

    ''' <summary>
    ''' Computes the R_D from Carlson symmetric form
    ''' </summary>
    ''' <param name="X"></param>
    ''' <param name="Y"></param>
    ''' <param name="Z"></param>
    ''' <returns>Construced from R_J(x,y,z,z) which is equal to R_D(x,y,z)</returns>
    ''' <remarks>http://en.wikipedia.org/wiki/Carlson_symmetric_form#Series_Expansion</remarks>
    Private Function RD(ByVal X As Double, ByVal Y As Double, ByVal Z As Double) As Double
        Dim sum, fac, lamda, dx, dy, dz, A, MinError As Double
        MinError = 0.000000001

        sum = 0
        fac = 1

        Do
            lamda = Math.Sqrt(X * Y) + Math.Sqrt(Y * Z) + Math.Sqrt(Z * X)
            sum += fac / (Math.Sqrt(Z) * (Z + lamda))

            fac /= 4

            X = (X + lamda) / 4
            Y = (Y + lamda) / 4
            Z = (Z + lamda) / 4

            A = (X + Y + 3 * Z) / 5

            dx = (1 - X / A)
            dy = (1 - Y / A)
            dz = (1 - Z / A)

        Loop While Math.Max(Math.Max(Math.Abs(dx), Math.Abs(dy)), Math.Abs(dz)) > MinError

        Dim E2, E3, E4, E5, result As Double
        E2 = dx * dy + dy * dz + 3 * dz ^ 2 + 2 * dz * dx + dx * dz + 2 * dy * dz
        E3 = dz ^ 3 + dx * dz ^ 2 + 3 * dx * dy * dz + 2 * dy * dz ^ 2 + dy * dz ^ 2 + 2 * dx * dz ^ 2
        E4 = dy * dz ^ 3 + dx * dz ^ 3 + dx * dy * dz ^ 2 + 2 * dx * dy * dz ^ 2
        E5 = dx * dy * dz ^ 3

        'http://dlmf.nist.gov/19.36#E2
        result = (1 - (3 / 14) * E2 + (1 / 6) * E3 + (9 / 88) * E2 ^ 2 - (3 / 22) * E4 - (9 / 52) * E2 * E3 + (3 / 26) * E5 - (1 / 16) * E2 ^ 3 + (3 / 40) * E3 ^ 2 + (3 / 20) * E2 * E4 + (45 / 272) * E2 ^ 2 * E3 - (9 / 68) * (E3 * E4 + E2 * E5))

        result = 3.0 * sum + fac * result / (A * Math.Sqrt(A))
        Return result

    End Function

Which make the implementation of the incomplete elliptic integrals of the first and second kind a piece of cake.

    ''' <summary>
    ''' Returns the imcomplete elliptic integral of the first kind 
    ''' </summary>
    ''' <param name="angle">In degrees, valid value range is from 0 to 90</param>
    ''' <param name="k">This function thakes k^2 as the parameter</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Function IncompleteEllipticIntegralFirstKind(ByVal angle As Double, ByVal k As Double) As Double
        Dim result As Double
        Dim ang As Double = Math.PI * angle / 180
        result = Math.Sin(ang) * RF(Math.Cos(ang) ^ 2, 1 - k * Math.Sin(ang) ^ 2, 1)
        Return result
    End Function

    ''' <summary>
    ''' Returns the imcomplete elliptic integral of the second kind 
    ''' </summary>
    ''' <param name="angle">In degrees, valid value range is from 0 to 90</param>
    ''' <param name="k">This function thakes k^2 as the parameter</param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Function IncompleteEllipticIntegralSecondKind(ByVal angle As Double, ByVal k As Double) As Double
        Dim result As Double
        Dim ang As Double = Math.PI * angle / 180
        result = Math.Sin(ang) * RF(Math.Cos(ang) ^ 2, 1 - k * Math.Sin(ang) ^ 2, 1) - (1 / 3) * k * Math.Sin(ang) ^ (3) * RD(Math.Cos(ang) ^ 2, 1 - k * Math.Sin(ang) ^ 2, 1)
        Return result
    End Function

I have tested the functions against the values given by Wolfram and they seem to match perfectly, However you cant do too much testing on a function as this, and I wont make any promises that it would give accurate results for any values. 

Elliptic curves    

The term Elliptic curves are in reallity misleading, as they have very little to do with ellipsis in practice. The name stems from Jacobi, who found the elliptic curves when inverting the elliptic integral of the first kind, called Jacobi's elliptic functions. Later Karl Weierstrass expanded the theory and found a simple elliptic function, expressed in a general form. 

They are quite difficult to explain so I am just going to mention the many implications of these curves. They are used by Wiles proof of Fermat's last theorem, and are also used in several problems in pure mathematics. Elliptic curves are also used in cryptography, as it is a very efficient encrypter, meaning it is often used in small devices as mobile phones etc.


References   

Online articles on elliptic functions, curves and equations:

Online pdf and power point documents  

Books   

Online books 

Videos   

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Kenneth Haugland
Engineer
Norway Norway
I hope that you like the stuff I have created and if you do wish to say thank you then a donation is always appreciated.
You can donate here[^].

Comments and Discussions

 
QuestionRamanujan type- Elliptical series1 PinmemberMember 1074252811-Apr-14 5:38 
AnswerRe: Ramanujan type- Elliptical series1 PinprofessionalKenneth Haugland2-Sep-14 4:14 
Questionangle bigger than pi/2 ? PinmemberMember 1012835226-Jun-13 18:25 
AnswerRe: angle bigger than pi/2 ? PinmemberMember 1012835226-Jun-13 19:19 
AnswerRe: angle bigger than pi/2 ? PinprofessionalKenneth Haugland26-Jun-13 21:29 
GeneralMy vote of 5 PinmemberMihai MOGA12-Apr-13 19:04 
GeneralRe: My vote of 5 PinprofessionalKenneth Haugland14-Apr-13 20:35 
GeneralMy vote of 5 Pinmemberhillsh_vi3-Apr-13 5:10 
GeneralRe: My vote of 5 PinmemberKenneth Haugland3-Apr-13 5:14 
GeneralMy vote of 5 PinmemberNik Steel2-Apr-13 10:42 
GeneralRe: My vote of 5 PinmemberKenneth Haugland2-Apr-13 19:18 
GeneralMy vote of 5 PinmemberPaul Conrad31-Mar-13 17:57 
GeneralRe: My vote of 5 PinmemberKenneth Haugland31-Mar-13 18:14 
GeneralRe: My vote of 5 PinmemberKenneth Haugland31-Mar-13 20:48 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web02 | 2.8.141022.2 | Last Updated 6 Apr 2013
Article Copyright 2013 by Kenneth Haugland
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid