This post has already been read 4044 times!

### 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:

Shape Cartesian 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):

### 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:

''' ''' Return the Complete Elliptic integral of the 1st kind ''' '''K(k^2) absolute value has to be below 1 ''' ''' Abramowitz and Stegun p.591, formula 17.3.11 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.

''' ''' Return the Complete Elliptic integral of the 2nd kind ''' '''E(k^2) absolute value has to be below 1 ''' ''' Abramowitz and Stegun p.591, formula 17.3.12 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:

''' ''' Computes the R_F from Carlson symmetric form ''' ''' ''' ''' ''' ''' http://en.wikipedia.org/wiki/Carlson_symmetric_form#Series_Expansion 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:

''' ''' Computes the R_D from Carlson symmetric form ''' ''' ''' ''' ''' Construced from R_J(x,y,z,z) which is equal to R_D(x,y,z) ''' http://en.wikipedia.org/wiki/Carlson_symmetric_form#Series_Expansion 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.

''' ''' Returns the imcomplete elliptic integral of the first kind ''' '''In degrees, valid value range is from 0 to 90 '''This function thakes k^2 as the parameter ''' ''' 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 ''' ''' Returns the imcomplete elliptic integral of the second kind ''' '''In degrees, valid value range is from 0 to 90 '''This function thakes k^2 as the parameter ''' ''' 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.

### Sources

Download source code for elliptic integrals in C# - 26.8 KB

Download source code and demo for elliptic integrals in C# - 47.8 KB

### References

Online articles on elliptic functions, curves and equations:

- http://mathworld.wolfram.com/EllipticIntegral.html
- http://mathworld.wolfram.com/EllipticCurve.html
- http://en.wikipedia.org/wiki/Torus
- http://en.wikipedia.org/wiki/Mordell%E2%80%93Weil_theorem
- http://www.math.hmc.edu/funfacts/ffiles/10006.3.shtml
- http://en.wikipedia.org/wiki/Ellipsoid
- http://en.wikipedia.org/wiki/Ellipse
- http://en.wikipedia.org/wiki/Modular_angle
- http://www.oocities.org/web_sketches/ellipse_notes/ellipse_arc_length/ellipse_arc_length.html
- http://www.math.wpi.edu/IQP/BVCalcHist/calc1.html
- http://en.wikipedia.org/wiki/Pendulum_%28mathematics%29#Arbitrary-amplitude_period
- Solving the simple pendulum by Jacobi elliptic function by Kenneth Shum