The Complete Elliptic Integrals and APL

by Joseph De Kerf

Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 102

In a previous note [1], a definition was given of the so-called common mean of two numbers, and it was shown how easy it is to program the algorithm in APL. In this note it is shown how the concept may be used to calculate the complete elliptic integrals of the first and second kind [2] [3]. User-defined functions in APL are given. Incomplete integrals will be treated in a following and final note.

The complete elliptic integrals of the first kind \(K(p)\) and of the second kind \(E(p)\) may be defined as:

\[\begin{align*} K(p) &= \int_0^{\pi/2} \frac{d\theta}{\sqrt{1-p^2\sin^2\theta}} \\[2ex] E(p) &= \int_0^{\pi/2} \sqrt{1-p^2\sin^2\theta}\;\; d\theta \end{align*}\]

with \(p^2 \le 1\) or \(-1 \le p \le 1\). Both complete elliptic integrals are even:

\[\begin{equation*} K(-p) = K(p) \text{ and } E(-p) = E(p) \end{equation*}\]

Complete elliptic integrals of the complementary modulus \(q\) are defined as \(K'(p) = K(q)\) and \(E'(p) = E(q)\), with \(q^2 = 1 - p^2\). The two auxiliary functions are also both even:

\[\begin{equation*} K'(-p) = K'(p) \text{ and } E'(-p) = E'(p) \end{equation*}\]

The complete elliptic integrals may be calculated for instance by numerical integration or by series expansion [2] [3]. Much easier to program, and this is especially true in APL, is the algorithm based on the common mean concept [2] [3]. Indeed, the complete elliptic integral of the first kind \(K(p)\) may be calculated by application of the relation:

\[\begin{equation*} K(p) = \frac{\pi}{2c(1,q)} \end{equation*}\]

where \(c(1,q)\) is the common mean of 1 and \(q = \sqrt{1-p^2}\)

A user-defined function CEI1, based on this procedure, may be:

      ∇CEI1[⎕]∇
    ∇R←CEI1 P
[1]   R←1,(1-P*2)*0.5
[2]   LAB:→(≠/R←(0.5×+/R),(×/R)*0.5)/LAB
[3]   R←(○1)÷+/R
    ∇

which for \(p = 0.0(0.1)0.9\) and comparison tolerance 1E¯15 gives:

      R←CEI1¨ P←0,0.1×⍳9
      3 1 15 10 20 15⍕P,R,[1.5] R
  0.0    1.5707963268    1.570796326794897
  0.1    1.5747455615    1.574745561517356
  0.2    1.5868678475    1.586867847454166
  0.3    1.6080486199    1.608048619930513
  0.4    1.6399998659    1.639999865864511
  0.5    1.6857503548    1.685750354812596
  0.6    1.7507538029    1.750753802915752
  0.7    1.8456939984    1.845693998374723
  0.8    1.9953027777    1.995302777664729
  0.9    2.2805491384    2.280549138422770

Results are displayed with 10 and 15 digits after the decimal point respectively. It should be noticed that for \(p = \pm 1\), the complementary modulus \(q\) is 0, and the complete elliptic integral of the first kind is infinite (\(\infty\)). A domain error is reported and the user-defined function is suspended. Eventually, the error report may be trapped and an appropriate reaction may be programmed.

The complete elliptic integral of the second kind \(E(p)\) may be calculated by application of the relation:

\[\begin{alignat*}{2} && E(p) &= \frac{K(p)}{2}\left( 2 - p^2 - \sum_{i=1}^{\infty} 2^i\bigl(x_i^2 - y_i^2\bigr) \right) \\[3ex] && E(p) &= \frac{K(p)}{2}\left( 2 - \sum_{i=0}^{\infty} 2^i\bigl(x_i^2 - y_i^2\bigr) \right) \\[3ex] &\text{or}\quad & E(p) &= \frac{\pi}{4c(1,q)}(2-S) \qquad\text{and}\qquad S = \sum_{i=0}^{\infty} 2^i\bigl(x_i^2 - y_i^2\bigr) \end{alignat*}\]

where \(c = (1,q)\) is the common mean of 1 and \(q = \sqrt{1-p^2}\), and the \(x_i\)s and \(y_i\)s are the successive arithmetic and geometric means while calculating this common mean:

\[\begin{equation*} \begin{aligned} & x_0 = 1 \\ & x_{i+1} = (x_i + y_i)/2 \end{aligned} \qquad\qquad\qquad\text{and}\qquad\qquad\qquad \begin{aligned} & y_0 = q \\ & y_{i+1} = \sqrt{x_i y_i} \end{aligned} \end{equation*}\]

A user-defined function CEI2, based on this procedure, may be:

      ∇CEI2[⎕]∇
    ∇R←CEI2 P;S;I
[1]   →(1=|P)/0,R←1
[2]   R←(1+S←1+I←¯1),(1-P*2)*0.5
[3]   LAB:S←S+(-/R*2)×2*I←I+1
[4]   →(≠/R←(0.5×+/R),(×/R)*0.5)/LAB
[5]   R←((○1)÷2×+/R)×2-S
    ∇

which for \(p = 0.0(0.1)1.0\) and comparison tolerance 1E¯15 gives:

      R←CEI2¨ P←0,0.1×⍳10
      3 1 15 10 20 15⍕P,R,[1.5] R
  0.0    1.5707963268    1.570796326794897
  0.1    1.5668619420    1.566861942021668
  0.2    1.5549685462    1.554968546242529
  0.3    1.5348334649    1.534833464923249
  0.4    1.5059416124    1.505941612360040
  0.5    1.4674622093    1.467462209339427
  0.6    1.4180833944    1.418083394448724
  0.7    1.3556611356    1.355661135571956
  0.8    1.2763499432    1.276349943169906
  0.9    1.1716970528    1.171697052781614
  1.0    1.0000000000    1.000000000000000

It should be noticed that for \(p = \pm 1\) the complementary modulus \(q\) is 0, the common mean \(c = (1,q)\) is 0, while the complete elliptic integral of the second kind is 1:\(E(\pm 1) = 1\). As such, in the first line of the function, for \(p = \pm 1\) the result \(R\) is set to 1 and flow of execution returns to the invoking expression. In the second line, the partial sum \(S\) and the counter \(I\) are initialised to 0 and -1 respectively. The loop, lines 3 and 4, may be written in one line of course, but the iteration has been split to increase readability, while the overhead is negligible. In the last line \(K(p)/2\) is multiplied by \(2-S\), giving finally \(R = E(p)\).

For both user-defined functions, the relative accuracy of the common means calculated is at least the current value of the comparison tolerance, for the results shown 1E¯15. In practice, however, this relative accuracy is slightly better, as the exact value of the common mean lies between the two iterations calculated and the result is set to the arithmetic mean of the final values of those iterations.

Note: The number of iterations needed and the cpu time depend on the current value of the comparison tolerance and on \(p\). Results of a benchmark for 1E¯15 and \(p = 0.1(0.2)0.9\) are shown below. The cpu times are in milliseconds.

\(p\)	0.1	0.3	0.5	0.7	0.9
Number of iterations	3	3	4	4	5
CEI1	112	112	137	137	163
CEI2	208	208	255	255	303

The cpu times are rather high, but are considerably better than for other methods such as for instance numerical integration and series expansion.

The calculations and the benchmark have been done on a MicroVAX 2000, with VAX APL Version 3.1, under VMS Version 4.7 [4]. The default value for the comparison tolerance in VAX APL is 1E¯15. The number of iterations has been measured by globalising the counter \(I\). The cpu times have been measured with the system function ⎕MONITOR.

References¶

[1] J. De Kerf: The Common Mean and APL. Vector Vol.11 No.3, pp21-23. 1995

[2] M. Abramowitz and I.A. Stegun (Eds): Handbook of Mathematical Functions. Applied Mathematics Series No. 55, National Bureau of Standards, Washington D. C., 1964 and subsequent revisions.

[3] J. Spanier and K.B. Oldham: An Atlas of Functions. Hemisphere Publishing Corporation, New York, New York 1987.

[4] VAX APL Reference Manual – Vols. 1 and 2 (Software Version V.3.0). Publication Numbers: AI-P142D-TE and AI-GV09B-TE. Digital Equipment Corporation, Maynard, Massachusetts, June 1987.

Production Manager’s Observation¶

When run in Dyalog/W, the sample evaluations produce very sightly different results! Both are run in the same comparison tolerance!

First Example …

      ⎕CT←1E¯15
      R←CEI1¨ P←0,0.1×⍳9
      3 1 15 10 20 15⍕P,R,[1.5] R
  0.0    1.5707963268    1.570796326794897
  0.1    1.5747455615    1.574745561517356
  0.2    1.5868678475    1.586867847454166
  0.3    1.6080486199    1.608048619930513
  0.4    1.6399998659    1.639999865864511
  0.5    1.6857503548    1.685750354812596
  0.6    1.7507538029    1.750753802915753
  0.7    1.8456939984    1.845693998374724
  0.8    1.9953027777    1.995302777664730
  0.9    2.2805491384    2.280549138422770

Second Example …

      ⎕CT←1E¯15
      R←CEI2¨ P←0,.1×⍳10
      3 1 15 10 20 15⍕P,R,[1.5] R
0.0    1.5707963268    1.570796326794897
0.1    1.5668619420    1.566861942021668
0.2    1.5549685462    1.554968546242529
0.3    1.5348334649    1.534833464923250
0.4    1.5059416124    1.505941612360040
0.5    1.4674622093    1.467462209339424
0.6    1.4180833944    1.418083394448724
0.7    1.3556611356    1.355661135571955
0.8    1.2763499432    1.276349943169906
0.9    1.1716970528    1.171697052781616
1.0    1.0000000000    1.000000000000000

Perhaps someone with Mathematica could say which is correct!!