15,038,345 members
Articles / Programming Languages / Visual Basic
Article
Posted 31 Mar 2013

55.7K views
26 bookmarked

# Elliptic integrals

Rate me:
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:

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

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:

VB.NET
''' <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.

VB.NET
''' <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:

VB.NET
''' <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:

VB.NET
''' <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.

VB.NET
''' <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

## Share

 Engineer Norway
No Biography provided

 First Prev Next
 partial arc length? scottfe6-May-15 14:29 scottfe 6-May-15 14:29
 Re: partial arc length? Kenneth Haugland6-May-15 20:08 Kenneth Haugland 6-May-15 20:08
 My vote of 5 phil.o27-Mar-15 10:05 phil.o 27-Mar-15 10:05
 Re: My vote of 5 Kenneth Haugland28-Mar-15 11:06 Kenneth Haugland 28-Mar-15 11:06
 Ramanujan type- Elliptical series1 Member 1074252811-Apr-14 5:38 Member 10742528 11-Apr-14 5:38
 Re: Ramanujan type- Elliptical series1 Kenneth Haugland2-Sep-14 4:14 Kenneth Haugland 2-Sep-14 4:14
 angle bigger than pi/2 ? Member 1012835226-Jun-13 18:25 Member 10128352 26-Jun-13 18:25
 Re: angle bigger than pi/2 ? Member 1012835226-Jun-13 19:19 Member 10128352 26-Jun-13 19:19
 Re: angle bigger than pi/2 ? Kenneth Haugland26-Jun-13 21:29 Kenneth Haugland 26-Jun-13 21:29
 My vote of 5 Ștefan-Mihai MOGA12-Apr-13 19:04 Ștefan-Mihai MOGA 12-Apr-13 19:04
 Re: My vote of 5 Kenneth Haugland14-Apr-13 20:35 Kenneth Haugland 14-Apr-13 20:35
 My vote of 5 hillsh_vi3-Apr-13 5:10 hillsh_vi 3-Apr-13 5:10
 Re: My vote of 5 Kenneth Haugland3-Apr-13 5:14 Kenneth Haugland 3-Apr-13 5:14
 My vote of 5 Nik Steel2-Apr-13 10:42 Nik Steel 2-Apr-13 10:42
 Re: My vote of 5 Kenneth Haugland2-Apr-13 19:18 Kenneth Haugland 2-Apr-13 19:18