The derivation of algorithms to compute elliptic integrals of the first and second kind by Landen-transformation

Desenvolvimento de algoritmos para o cálculo de integrais elípticas de primeira e segunda ordens por meio da Transformação de Landen

ARTIGOS

The derivation of algorithms to compute elliptic integrals of the first and second kind by Landen-transformation

Desenvolvimento de algoritmos para o cálculo de integrais elípticas de primeira e segunda ordens por meio da Transformação de Landen

Norbert Rösch

Geodätisches Institut, Universität Karlsruhe, Englerstrasse 7, D-76128 Karlsruhe, Germany, [e-mail]

ABSTRACT

This paper deals with elliptic integrals of first and second kind and their solution by Landen-transformation. It is explained how the Landen-transformation works and how the different algorithms can be used in geodetic problem solving. The algorithms are recursive and thus easy to implement and stable. Further on the testing is quite simple. The numerical results of the Landen-transformation are compared to this computed with Maple V to show the exactness of the algorithms. For two algorithms the source-code is depicted.

Keywords: Arc of Meridian; Numerical Integration; Elliptic Integral; Landen-Transformation.

RESUMO

Neste trabalho enfoca-se as integrais elípticas de primeira e segunda ordens e sua solução por meio da Transformação de Landen. Explica-se como a Transformação de Landen opera e como diferentes algoritmos podem ser usados na solução de problemas em Geodésia. Os algoritmos são recursivos e por isso, fáceis de implementar. Testar os algoritmos também é fácil. Os resultados da Transformação de Landen são comparados aos resultados obtidos com Maple V para mostrar a exatidão dos algoritmos. Apresenta-se o código fonte para os dois algoritmos.

Palavras-Chave: Arco de Meridiano; Integração Numérica; Integral Elíptica; Transformação de Landen.

1. INTRODUCTION

In 1771 and 1775 the British mathematician John Landen (1719-1790) published an appropriate (which means applicable at this time) method to solve elliptic integrals (see LANDEN (1771), LANDEN (1775)). This approach is now called the Landen-transformation or Landen-substitution respectively. Remarkably the above mentioned transformation was published before Joseph Liouville (1809-1882) proved that there exists no anti derivative for the elliptic integrals. However, the definite integral can be solved.

In addition to the Landen-transformation there are still other approaches to solve the elliptic integrals. These are for example the series expansion or the numerical integration. Both topics comprise a number of methods. Also in the group of the transformation/substitution methods another kind of transformation is known, the so-called Gauss-transformation. However, these approaches are beyond the scope of this paper.

In fact there are different ways to depict the elliptic integrals but most commonly the notation of Legendre is used. According to this notation the definition of the incomplete elliptic integrals is:

Elliptic integral of first kind

Elliptic integral of second kind

The terms k and Φ are denoted as modulus and amplitude respectively.

In mathematical and cartographic geodesy there are different problems where these integrals are found. Only a few examples of these should be shown in this introduction. First the computation of the elliptic arc length S should be mentioned which leads to (see HECK, (1995))

or

(a = semi-major axis of the ellipsoid of revolution; e = first eccentricity, Φ₁ = latitude) where the first term of the formula includes an elliptic integral of second kind with the modulus e. Hopfner explained in (HOPFNER,(1938)) how the arc length can be computed by different methods. In this context, he also mentioned the Landen-transformation.

Especially formula 3 reveals that the integral has to be computed with a high accuracy, because the result of E(e, Φ₁) is multiplied with the semi-major axis of the ellipsoid of revolution. So the error of the computation of the elliptic integral should be smaller than 1.10^-11 to get the arc length with an accuracy of 1 mm. This is the reason why for all further computations a threshold ε for the recursive processes of ε < 1.10^-15 is introduced.

Also the transformation between isometric and geographic coordinates could be computed with elliptic integrals of the second kind (HECK, 1995, p. 262). Thus, their solution is still of interest and the transformation of isometric coordinates is their main area of application.

Nevertheless there are other geodetic problems which lead to elliptic integrals. For example the so-called indirect geodetic problem. In this context the length of the geodesic s can be calculated by:

2. THE COMPUTATION OF ELLIPTIC INTEGRALS OF THE FIRST KIND

First of all have a look at the definite elliptic integral of first kind given by:

The substitution of Φ by ω leads to ¹ :

and the derivate:

and furthermore we obtain:

The division of (7) and (8) combined with the substitution of cos(2ω) leads to

Now the above shown elliptic integral of first kind can be written as:

Both sides of the formula are elliptic integrals of the first kind. But whereas the left integral is original the right integral is of different modulus and amplitude. This is abbreviated by:

Note: Due to κ < 1 we obtain: (1 + κ)² < 4 or . As.

From formula (7) we get:

sin Φ(κ + cos 2ω) = sin 2ωcosΦ

or

which yields:

Thus we got an easy way to use formula to compute ω. With sin(2ω - Φ) < sin Φ we obtain ω < Φ. I.e. in equation (11) the modulus of the left integral is smaller whereas the amplitude is bigger. This can be repeated and as a result increasing module and decreasing amplitudes are obtained. Written as a chain of elliptic integrals of the first kind with increasing module we get (κ = κ₀) and (Φ = Φ₀):

with:

and referring to equation (13):

Due to:

equation (14) could be written as:

With this equation we gained a practical formula to compute the values of an elliptical integral of first kind.

But there is still another way to read equation (11). It follows directly:

which means that with a given q and ω an elliptical integral with a smaller modulus is found on the right. The relations between p, Φ and q, ω are:

and

Written as a suitable formula the above could be combined to (with κ₀ = q and Φ₀ = ω)

It must be considered that

and thus

and referring to (22) we get

2.1 Solution with Increasing Modulus

In this chapter numerical aspects of an algorithm based on (18) are outlined. This formula could be easily implemented as recursive algorithm. At first, the convergence of this method should be examined. Therefore the module of each iteration loop is considered. The results are depicted in table 1 (the values are computed with Φ = 40º). Table (1) shows, that even for k₀ = 0.1 the algorithm needs only 5 steps to achieve the break criterion |κ_n-1 – κ_n| = ε < 1.10^-15 (the results in all the tables are calculated with ε = 1. 10^-15). This gives an idea for the convergence of the Landen-transformation.

Furthermore it is remarkable that for big module the break criterion of the iteration is fulfilled faster (as expected). For κ = 0.1 e. g. there are 5 steps needed, whereas for κ = 0.9 4 steps are enough. If the break criterion ε would be changed to ε = 1. 10 ^-14, than for κ = 0.9 even 3 steps instead of 4 are sufficient.

The following table shows the relation between the modulus and its correlated amplitude. As an example the modulus κ₀ = 0.1 is used (table 2). Whereas the modulus in column 2 increases from 0.1 to 1.0 within 5 steps, the amplitude (see column 3) decreases from 40.0 degrees to about 17 degrees. Table 3 gives an overview about the numerical results calculated by increasing module. These values could be compared with those of table 6. The implementation of the algorithm based on (18) is shown below. The source code is written in C. Like all recursive processes the function to compute F (k, Φ) is rather compact. Further on the example is intended to show how easy equation (18) can be implemented.

double F_inc(double modul,double amplitude, double eps){

amplitude = (asin(modul*sin(amplitude))+amplitude)/2.;

if ((1.-modul) < eps)

return log(tan(amplitude/2.+EllInt_PI/4.));

else

return 2./(1.+modul)*F_inc((2.*sqrt(modul))/(1.+modul),amplitude,eps); }

The quantities modul and amplitude in the call of the function are self-explanatory. eps denotes the break criterion for the recrusive process and can be any arbitrary value. By changing this quantity the user can control the processing time and the accuracy of the algorithm. The quantity EllInt_PI is a symbolic constant and simply stands for the circle constant π.

2.2 Solution with Decreasing Modulus

In this chapter an iteration based on (25) is outlined. As already shown before, this formula could be used to compute module converging to 0. To calculate the different values for the module and the amplitudes the formulas (20) and (21) could be applied.

According to what we said before, in this chapter the results of the tables 1 - 3 should now be computed with an algorithm, which works with decreasing module.

As column two in table 4 points out, this algorithm needs for κ = 0.1 only four steps to reach its break criterion, whereas for κ = 0.5 and κ = 0.9 five steps are needed. In chapter 2.1 this was contrariwise. The method using increasing module needed five steps for κ = 0.1 and for κ = 0.5 and κ = 0.9 four are sufficient. This result is exactly what is expected. Because the goal of the algorithm working with increasing module is to reach a modulus of 1, whereas the other one should reach a modulus of 0.

Now have a look at table 5, which depicts the decreasing module and the increasing amplitudes. In the last row of column three the angle is about two times bigger than that one of the last but one iteration. Besides the angle of this iteration is almost twice as big as the angle computed in the step before. Of cause this is not incidental. The connection between Φ_n and Φ_n+1 is given by

and thus

Hence, for small module the angle will double from one iteration to the next.

Finally the results of table 3 and table 6 should be compared. In seven out of twenty cases the results are identical. Nine times the difference is about 1.10^-15. In three cases the difference is 2.10^-15 and only once it is 11.10^-15. This means that the results confirm each other.

2.3 Solution with Decreasing Modulus and Series Expansion

Using decreasing module there is still another way to find an algorithm based on formula (22). As table 4 shows, the module decrease rapidly. Even for κ = 0.9 the modulus of step three is less than 5.10^-4. This seems to be small enough to try a series expansion.

The term (1 – x)^-1/2 (|x| < 1 but it works better if |x| << 1 could be written as

Applied on F (k, Φ) we obtain

Thus (22) will become to

This formula could be also implemented as a recurrence. The break criterion would then be a modulus k which is smaller than 1.10^-3. An approximation leads us to the following results. If κ = 1.10^-3 (worst case scenario) we get for the 4^th term R in equation (30).

The effect of this term is marginal. So it is sufficient to take the first three terms of the series in formula (30). Table 7 depicts the results computed with the first three terms.

2.4 Comparison of the Results

Obviously it is difficult to decide which of the results is right. Therefore all computations are compared with those gained by Maple V, which are treated as a reference. The results are illustrated in table 8. Comparing table 8 with the tables 3 - 7 shows that only two out of twenty values have a difference of more than 1.10^-15. One of these differences is 2.10^-15 and the other is 8.10^-15, computed for F (0.999, 90º) with increasing module.

   

At first sight, the deviation d = 8.10 ^-15 seems to be a contradiction to the kind of method which is applied. Because increasing module should lead to rather good results, as the number of iterations should be small. If we look at the module, which are κ₀ = 0.999, κ₁ = 0.999999874874898, κ₂ = 0.999999999999998 and κ₃ = 1.0, we realize, that

which means that the break criterion is failed only by 1.10^-15. I. e. even the last but one iteration κ₂ almost fulfills the break criterion. Due to numerical aspects the next iteration delivers no further contribution to the result. So the derivation between the reference and the computed solution is relatively big. Nevertheless the Landen-transformation is a useful method to compute elliptic integrals of the first kind. Even the biggest difference of 8.10^-15 between the two computed results easily fulfills the geodetic requirements mentioned in section 1.

3. THE COMPUTATION OF ELLIPTIC INTEGRALS OF SECOND KIND

Also for the elliptic integral of the second kind (see 2) there exist explicit solutions for two specific modules. The first is κ = 1 and the second is κ = 0.

With this, the same substitution as in section 2 is performed. From (6) we derive

and thus

and further with (7) combined with (34)

With the analogous conversions as in (9) we get

Some more derivations of the right side are leading to

and

Note:

From this we obtain

Changing the indexes to κ = κ₀ and q = κ₁ as well as Φ = Φ₀ and ω = Φ₁ we obtain

In this formula we have on the left the original elliptic integral of the second kind with its modulus κ₀ and its amplitude Φ₀, whereas on the right we obtain an elliptic integral of the second kind with a smaller modulus κ₁ combined with an elliptic integral of first kind whereas the amplitude is modified in both cases to Φ₁.

3.1 Solution with Increasing Modulus

In fact there are two different ways to read formula (42). The first one is that by subsequent conversion of the term E(κ_n, Φ_n) until

and thus

is obtained. Obviously this can be implemented as a recursive algorithm too. It must be considered that this formula encompasses two recursive processes because both, the elliptic integral of second kind and that of the first kind have to be calculated that way. But for the computation of the latter we already found various solutions in section 2.

Further on we can apply (19) together with (20) and find

In this formula the elliptic integral of the first kind has the same modulus and the same amplitude as the elliptic integral of the second kind. This conversion could be repeated in the same way as we have done before. As a result we get

The above shown equation seems to be not so time consuming than (46) because the elliptic integral of the first kind has to be computed only once. Due to numerical aspects it is not possible to implement it as a recursive algorithm. In step n the modulus is approximately 1, because this is the break criterion. But now, beginning with this modulus, the decreasing chain of the moduli has to be computed. This of cause is also a recursive process (at least in theory). In fact we get as a result of

As a consequence the recursive process will have no end. This problem could only be solved with a quasi – recursive algorithm, which puts the computed moduli of step i in a stack from which the product of the moduli is derived (see mark '!' in equation (46)). If done so, the numerical values are exactly the same as the results computed with (44) (see table 9 together with table 10).

As the results are nearly the same, the algorithm based non (46) is about 35% faster than this based on (44)² . However the implementation of (46) is much more complicated, because it is not strictly recursive, as already mentioned.

An implementation of an algorithm based on (44) is depicted below.

double E_inc(double modul,double amplitude, double eps) {

double n_amp,n_mod;

n_amp = (asin(modul*sin(amplitude))+amplitude)/2.;

if ((1.-modul)<eps)

return (1+modul)*sin(n_amp)+(1-modul)*

F_inc((2.*sqrt(modul))/(1.+modul),n_amp,eps)-modul*sin(amplitude);

else {

n_mod = (2.*sqrt(modul))/(1.+modul);

return (1.+modul)*E_inc(n_mod,n_amp,eps)+(1.-modul)*

F_inc(n_mod,n_amp,eps)-modul*sin(amplitude); } }

The example shows further that the recursive process invoked by the call of E_inc³ comprises a second one invoked by F_inc (see section 2.1). I. e. the second is nested within the first. Thus the algorithm is extremely compact.

3.2 Solution with Decreasing Modulus

An algorithm to compute E(p, Φ) with decreasing module based on (42) is obtained by

In this formula we assume that E(κ₁, Φ₁) has to be computed and E(κ₀, Φ₀) is known. Therefore it makes sense to change the indexes and we get

The terms κ₁ and Φ₁ are computed with (20) and (21) respectively. So the modulus κ₁ is smaller than κ₀ whereas Φ₁ is larger than Φ₀. E(κ_1, Φ₁) could be treated in the same way as E(κ_0, Φ₀) until we get κ < ε, which means in practice κ ≈ 0. With

a break criterion is given. Thus a recursive process to compute E(p, Φ) is found. The numerical results are shown in table 11.

In (49) the term (1 – κ₁)F(κ_1, Φ₁) could be replaced in the same way as we have done in section 3.1. Now the elliptic integral of the first kind needs to be computed only once.

3.3 Solution with Decreasing Modulus and Series Expansion

In the same way we did in section 2.3, we can combine the Landen-transformation with the series expansion. In this context this seems to be even more advantageous because the computation of an elliptic integral of the second kind is even more time-consuming.

The series expansion of

yields

The substitution of x = κ²sin²Φ applied to (2) leads to

As above mentioned formula (51) will deliver the best results if |x| << 1and thus κ << 1. Hence we decrease the modulus by Landen-transformation to apply the series expansion to the transformed integral subsequently.

If we look at the 4^th term of the series in formula (53) and integrate it assuming k = 0.001 and the maximum amplitude, the result r is r < 2.7.10^-26. So it is sufficient to take only the first three terms of the series. Thus the fourth and all further terms are omitted. So formula (49) in combination with (53) yields

This formula only depicts the last term of the recursion. At this stage the modulus κ is κ < 10^-3.

It makes sense to compute the elliptic integral of first kind in equation (54) also by decreasing module in combination with the series expansion, because in each step of the recursion the modulus decreases and thus the break criterion is reached faster.

3.4 Comparison of theRresults

The comparison between the tables 9 - 12 shows that most of the values are identical with those of table 13, which is treated as a reference. The biggest difference concerns E(0.999, 90º), which leads to 2.10^-15 between table 12 and 13.

Beside the accuracy there is still another criterion which should not be neglected. This criterion is the processing time. All algorithms of section 3 comprise two recursive processes. So it should be considered to take the right combination of both, to get the best result. Whereas in this case best means fastest. E. g. it seems to be different if we solve E(κ, Φ) by decreasing module, then it makes sense to apply decreasing module too, in oder to compute the appropriate elliptic integral of the first kind.

This can easily be proved using E(0.5, 90º) as an example. The algorithm based on 54 in combination with that based on 30 is about 25% faster than that based on 54 together with 25. The same test applied on E(0.1, 90º) delivers still an advantage of 20% between the two approaches. So the decision between the different methods is mainly influenced by the processing time, because the accuracy is the same.

4 CONCLUSIONS

As an example the result of the arc length for Φ = 48º is computed. With (3) we get

S = 5317885.233 m

(Φ denotes the amplitude; further on the parameters of the ellipsoid of Bessel are chosen) which is exactly the same result as published by (SCHÖDLBAUER, 1981). The result remains the same regardless which algorithm is applied. Even the combination of decreasing and increasing module (e.g. formula 49 together with 18 or 42 together with 25) leads to the same result. Therefore it is almost needless to explain that as a consequence also the transformation from geographical coordinates to conformal coordinates (e. g. of type Gauss-Krueger or UTM) yields the same results.

The advantages of the algorithms based on the Landen-transformation are obvious. They are flexible which means that they can be applied to different problems and the accuracy can be fitted to the requirements. Further on they are simple to implement as a recursive process and thus the test of the implementation is easy as well. As a consequence the presented solutions are an easy to use alternative to the known approaches.

ACKNOWLEDGMENT

The author is very thankful to Dr. Michael Kern, he supported the first steps of this work. My very special thanks go to Mr. Reginald Atkins, who encouraged me to write this paper.

(Invited Paper. Recebido em fevereiro de 2011).

• HECK, B. Rechenverfahren und Auswertemodelle in der Landesvermessung. Number 2. Herbert Wichmann Verlag, Karlsruhe, 1995.
• HOPFNER, F. GPS Zur Berechnung des Meridianbogens, Z. f. Verm.wesen, nş, 1938, p. 620627.
• LANDEN, J. A disquisition concerning certain fluents, which are assignable by the arcs of the conic sections; wherein are investigation some new and useful theorems for computing such fluents. Philosophical TransactionsRoyal Society of the Royal Society of London, 1771, p. 298309.
• LANDEN, J. An investigation of a general theorem for finding the length of any arc of any conic hyperbola, by means of two elliptic arcs, with some other new and useful theorems deduced therefrom. Philosophical TransactionsRoyal Society of the Royal Society of London, 1775, p. 283 289.
• SCHLÖMILCH, O. Compendium der Höheren Analysis Braunschweig, 1866.
• SCHÖDLBAUER, A. Rechenformeln und Rechenbeispiele zur Landesvermessung, volume 1. Herbert Wichmann Verlag, Karlsruhe, 1981.

1

In this section as well as in section 3 all explanations concerning the derivation of the Landen-transformation are based on (SCHLÖMILCH, 1866).

2

The processing time depends mainly on the modulus of cause. So its difficult to give a precise estimation.

3

The meaning of the different quantities is the same as already described in section 2.1 for F_inc

Publication Dates

• Publication in this collection
  02 Sept 2011
• Date of issue
  Mar 2011

History

• Received
  Feb 2011