# Complex incomplete elliptic integral of the second kind

Posted on January 30, 2020 by Stéphane Laurent

The incomplete elliptic integral of the second kind is usually defined by $E(\phi \mid k) = \int_0^\phi \sqrt{1 - k^2\sin^2t}\,\mathrm{d}t$ for $$-\frac{\pi}{2} \leqslant \phi \leqslant \frac{\pi}{2}$$ and $$k^2 \leqslant \frac{1}{\sin^2\phi}$$.

It can be evaluated in R with the gsl package:

f <- function(t, k) sqrt(1 - k^2*sin(t)^2)
phi <- 1; k <- 0.5
integrate(f, 0, phi, k = k)$value ## [1] 0.9648765 gsl::ellint_E(phi, k) ## [1] 0.9648765 That said, this elliptic integral is more generally defined by $E(\phi \mid m) = \int_0^\phi \sqrt{1 - m\sin^2t}\,\mathrm{d}t$ for $$-\frac{\pi}{2} \leqslant \phi \leqslant \frac{\pi}{2}$$ and $$-\infty < m \leqslant \frac{1}{\sin^2\phi}$$. # Evaluation of $$E(\phi \mid m)$$ Thus, gsl allows to evaluate $$E(\phi \mid m)$$ for a positive value of the parameter $$m$$ only, with the ellint_E function. However it is possible to use gsl to evaluate $$E(\phi \mid m)$$ for a negative parameter $$m$$, with the help of the Carlson elliptic integrals $$R_F$$ and $$R_D$$. These elliptic integrals are evaluated in gsl by the functions ellint_RF and ellint_RD. E <- function(phi, m){ if(phi == 0){ 0 }else if(phi == pi/2 && m == 1){ 1 }else if(phi >= -pi/2 && phi <= pi/2){ sine <- sin(phi) sine2 <- sine*sine cosine2 <- 1 - sine2 oneminusmsine2 <- 1 - m*sine2 sine * (gsl::ellint_RF(cosine2, oneminusmsine2, 1) - m * sine2 * gsl::ellint_RD(cosine2, oneminusmsine2, 1) / 3) }else{ NaN } } Let’s try it: E(phi = 1, m = -5) ## [1] 1.493736 You can compare with EllipticE[1,-5] on WolframAlpha. Here is how the function looks like for a fixed value of $$\phi$$: opar <- par(mar = c(5,4,1,1)) phi <- pi/4 curve( E(phi, x), from = -5, to = 1/sin(phi)^2, lwd = 2, xlab = expression(italic(m)), ylab = expression(paste(italic(E), "(", italic(phi1), " | ", italic(m), ")")) ) par(opar) # Extending $$E(\phi \mid m)$$ to arbitrary real $$\phi$$ There is a continuous extension of $$E(\phi \mid m)$$. It is given by the formula $E(\phi + k\pi \mid m) = 2 k E\Bigl(\frac{\pi}{2} \mid m\Bigr) + E(\phi \mid m)$ for any $$k \in \mathbb{Z}$$. E <- function(phi, m){ if(phi == 0){ 0 }else if(phi == pi/2 && m == 1){ 1 }else if(phi >= -pi/2 && phi <= pi/2){ sine <- sin(phi) sine2 <- sine*sine cosine2 <- 1 - sine2 oneminusmsine2 <- 1 - m*sine2 sine * (gsl::ellint_RF(cosine2, oneminusmsine2, 1) - m * sine2 * gsl::ellint_RD(cosine2, oneminusmsine2, 1) / 3) }else if(phi > pi/2){ k <- 0 while(phi > pi/2){ phi <- phi - pi k <- k + 1 } 2*k*E(pi/2, m) + E(phi, m) }else{ k <- 0 while(phi < -pi/2){ phi <- phi + pi k <- k - 1 } 2*k*E(pi/2, m) + E(phi, m) } } Let’s try it: E(phi = 7, m = -5) ## [1] 12.25899 and compare with EllipticE[7,-5] on WolframAlpha. # Extending $$E(\phi \mid m)$$ to complex $$\phi$$ and $$m$$ Now we will construct the analytical continuation of $$E(\phi \mid m)$$. I found a nice implementation of the Carlson elliptic integrals here, which can easily be generalized to complex values of the parameters: RF <- function(x, y, z, minerror = 1e-7){ x <- as.complex(x); y <- as.complex(y); z <- as.complex(z) dx <- dy <- dz <- Inf while(max(dx,dy,dz) > minerror){ lambda <- sqrt(x)*sqrt(y) + sqrt(y)*sqrt(z) + sqrt(z)*sqrt(x) x <- (x + lambda) / 4 y <- (y + lambda) / 4 z <- (z + lambda) / 4 A <- (x+y+z) / 3 dx <- Mod(1 - x/A) dy <- Mod(1 - y/A) dz <- Mod(1 - z/A) } E2 <- dx*dy + dy*dz + dz*dx E3 <- dy*dx*dz (1 - E2/10 + E3/14 + E2*E2/24 - 3*E2*E3/44 - 5*E2*E2*E2/208 + 3*E3*E3/104 + E2*E2*E3/16) / sqrt(A) } RD <- function(x, y, z, minerror = 1e-7){ x <- as.complex(x); y <- as.complex(y); z <- as.complex(z) dx <- dy <- dz <- Inf s <- 0; fac <- 1 while(max(dx,dy,dz) > minerror){ lambda <- sqrt(x)*sqrt(y) + sqrt(y)*sqrt(z) + sqrt(z)*sqrt(x) s <- s + fac / (sqrt(z) * (z + lambda)) fac <- fac/4 x <- (x + lambda) / 4 y <- (y + lambda) / 4 z <- (z + lambda) / 4 A <- (x + y + 3*z) / 5 dx <- Mod(1 - x/A) dy <- Mod(1 - y/A) dz <- Mod(1 - z/A) } E2 <- dx*dy + dy*dz + 3*dz*dz + 2*dz*dx + dx*dz + 2*dy*dz E3 <- dz*dz*dz + dx*dz*dz + 3*dx*dy*dz + 2*dy*dz*dz + dy*dz*dz + 2*dx*dz*dz E4 <- dy*dz*dz*dz + dx*dz*dz*dz + dx*dy*dz*dz + 2*dx*dy*dz*dz E5 <- dx*dy*dz*dz*dz 3*s + fac * (1 - 3*E2/14 + E3/6 + 9*E2*E2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26 - E2*E2*E2/16 + 3*E3*E3/40 + 3*E2*E4/20 + 45*E2*E2*E3/272 - 9*(E3*E4 + E2*E5)/68) / A / sqrt(A) } An important note. In the link I gave above, the author takes $\lambda = \sqrt{xy} + \sqrt{yz} + \sqrt{zx}$ whereas I take $\lambda = \sqrt{x}\sqrt{y} + \sqrt{y}\sqrt{z} + \sqrt{z}\sqrt{x}.$ This is important: the square root of a product is not equal to the product of the square roots in general, for complex numbers! Now here is the implementation of the incomplete elliptic integral of the second kind allowing complex values of $$\phi$$ and $$m$$: E <- function(phi, m){ if(phi == 0){ 0 }else if(phi == pi/2 && m == 1){ 1 }else if(Re(phi) >= -pi/2 && Re(phi) <= pi/2){ sine <- sin(phi) sine2 <- sine*sine cosine2 <- 1 - sine2 oneminusmsine2 <- 1 - m*sine2 sine * (RF(cosine2, oneminusmsine2, 1) - m * sine2 * RD(cosine2, oneminusmsine2, 1) / 3) }else if(Re(phi) > pi/2){ k <- 0 while(Re(phi) > pi/2){ phi <- phi - pi k <- k + 1 } 2*k*E(pi/2, m) + E(phi, m) }else{ k <- 0 while(Re(phi) < -pi/2){ phi <- phi + pi k <- k - 1 } 2*k*E(pi/2, m) + E(phi, m) } } Let’s try it: E(phi = 7+2i, m = -5) ## [1] 7.763389+5.631026i and compare with EllipticE[7+2I,-5] on WolframAlpha. I also find the same results as the ones of WolframAlpha when $$\phi$$ and $$m$$ are complex. E(phi = 1+1i, m = -5i) # same as Wolfram ## [1] -0.66421+1.746623i E(phi = 7+2i, m = -5i) # same as Wolfram ## [1] 14.95287+4.83345i ## Resorting to the Appell function As a sanity check, let’s propose another way. This way consists in using the relation between the incomplete elliptic integral of the second kind and the first Appell hypergeometric function: $E(\phi, m) = \sin\phi \times F_1\Bigl(\frac12, \frac12, -\frac{3}{2}, \sin^2\phi, m\sin^2\phi\Bigr),$ for $$-\frac{\pi}{2} \leqslant \Re(\phi) \leqslant \frac{\pi}{2}$$. Unfortunately there’s no available R implementation of $$F_1$$ allowing complex values of the variables. Then we will use the Euler integral representation $F_1\Bigl(\frac12, \frac12, -\frac{3}{2}, z_1, z_2\Bigr) = \frac12\int_0^1 \frac{\sqrt{1-z_2t}}{\sqrt{t}\sqrt{1-z_1t}}\,\mathrm{d}t.$ integrand <- function(t, z1, z2){ sqrt(1-t*z2) / sqrt(t) / sqrt(1-t*z1) } integrandRe <- function(t, z1, z2) Re(integrand(t, z1, z2)) integrandIm <- function(t, z1, z2) Im(integrand(t, z1, z2)) E <- function(phi, m){ if(phi == 0){ 0 }else if(phi == pi/2 && m == 1){ 1 }else if(Re(phi) >= -pi/2 && Re(phi) <= pi/2){ sine <- sin(phi) sine2 <- sine*sine re <- integrate(integrandRe, 0, 1, z1 = sine2, z2 = m*sine2) im <- integrate(integrandIm, 0, 1, z1 = sine2, z2 = m*sine2) sine/2 * (re$value + 1i*im\$value)
}else if(Re(phi) > pi/2){
k <- 0
while(Re(phi) > pi/2){
phi <- phi - pi
k <- k + 1
}
2*k*E(pi/2, m) + E(phi, m)
}else{
k <- 0
while(Re(phi) < -pi/2){
phi <- phi + pi
k <- k - 1
}
2*k*E(pi/2, m) + E(phi, m)
}
}
#
E(1+1i, -5i) # almost the same as before
## [1] -0.664216+1.746635i
#
E(7+2i, -5i) # almost the same as before
## [1] 14.95287+4.83346i

The slight differences are due to the numerical integrations. The previous way is preferable.