--- title: "Complex incomplete elliptic integral of the second kind" author: "Stéphane Laurent" date: '2020-01-30' tags: R, maths, special-functions output: md_document: variant: markdown preserve_yaml: true html_document: highlight: kate keep_md: no prettify: yes linenums: yes prettifycss: minimal highlighter: pandoc-solarized editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, fig.path = "figures/EllipticII-") ``` 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: ```{r} f <- function(t, k) sqrt(1 - k^2*sin(t)^2) phi <- 1; k <- 0.5 integrate(f, 0, phi, k = k)$value gsl::ellint_E(phi, k) ``` 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`. ```{r} 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: ```{r} E(phi = 1, m = -5) ``` You can compare with [`EllipticE[1,-5]`](https://www.wolframalpha.com/input/?i=EllipticE%5B1%2C+-5%5D) on WolframAlpha. Here is how the function looks like for a fixed value of $\phi$: ```{r curveE, fig.width=5, fig.height=5} 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}$. ```{r} 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: ```{r} E(phi = 7, m = -5) ``` and compare with [`EllipticE[7,-5]`](https://www.wolframalpha.com/input/?i=EllipticE%5B7%2C+-5%5D) 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](https://www.codeproject.com/Articles/566614/Elliptic-integrals), which can easily be generalized to complex values of the parameters: ```{r} 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$: ```{r} 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: ```{r} E(phi = 7+2i, m = -5) ``` and compare with [`EllipticE[7+2I,-5]`](https://www.wolframalpha.com/input/?i=EllipticE%5B7%2B2I%2C+-5%5D) on WolframAlpha. I also find the same results as the ones of WolframAlpha when $\phi$ and $m$ are complex. ```{r} E(phi = 1+1i, m = -5i) # same as Wolfram E(phi = 7+2i, m = -5i) # same as Wolfram ``` ## 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. $$ ```{r} 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 # E(7+2i, -5i) # almost the same as before ``` The slight differences are due to the numerical integrations. The previous way is preferable.