Simultaneously computing a complete elliptic integral and its complement - MathOverflow most recent 30 from http://mathoverflow.net 2013-05-21T06:39:53Z http://mathoverflow.net/feeds/question/34363 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/34363/simultaneously-computing-a-complete-elliptic-integral-and-its-complement Simultaneously computing a complete elliptic integral and its complement J. M. 2010-08-03T09:42:43Z 2012-04-30T11:17:39Z <p>The complete elliptic integral of the first kind</p> <p><code>$K(m)=\int_0^{\pi/2}\frac{\mathrm{d}t}{\sqrt{1-m\sin^2t}}$</code></p> <p>is easily computed via the arithmetic-geometric mean iteration; to wit,</p> <p><code>$K(m)=\frac{\pi}{2M(1,\sqrt{1-m})}$</code></p> <p>where <code>$M(a,b)$</code> is the arithmetic-geometric mean of <code>$a$</code> and <code>$b$</code>. With a little more trickery, the iteration can be hijacked to compute the complete elliptic integral of the second kind <code>$E(m)$</code> as well.</p> <p>In a number of applications, it happens that one needs both the values of <code>$K(m)$</code> and its complement <code>$K(1-m)$</code> (and sometimes similarly for <code>$E(m)$</code> and <code>$E(1-m)$</code>).</p> <p>My question is, apart from having to do an AGM iteration for each of <code>$K(m)$</code> and <code>$K(1-m)$</code>, is there an algorithm (maybe a modification of the basic AGM iteration) that simultaneously generates both <code>$K(m)$</code> and its complement? I would also be interested in seeing also an extension of this algorithm, if one exists, for computing <code>$E(m)$</code> as well (after which <code>$E(1-m)$</code> is easily computed via Legendre's relation).</p> http://mathoverflow.net/questions/34363/simultaneously-computing-a-complete-elliptic-integral-and-its-complement/34457#34457 Answer by SandeepJ for Simultaneously computing a complete elliptic integral and its complement SandeepJ 2010-08-04T00:06:25Z 2010-08-04T14:28:24Z <p>There are two possible ways to attack this problem</p> <ol> <li><p>Both K and K' can be expressed in terms of the Theta function as described here <a href="http://mathworld.wolfram.com/EllipticModulus.html" rel="nofollow">http://mathworld.wolfram.com/EllipticModulus.html</a>. If you compute $\Theta_3$, you can get both at the same time.</p></li> <li><p>The other way is to observe that both K and K' are expressible in terms of the hypergeometric function $_2F_1(\frac{1}{2}, \frac{1}{2} ; 1; m)$. They are solutions of the same self-adjoint Gauss hypergeometric differential equation (since the equation is invariant under the transformation (m $\rightarrow$ 1-m))</p> <p>$(k^3 - k)\frac{d^2y}{dk^2} + (3k^2 -1)\frac{dy}{dk} + ky = 0 $ </p> <p>By virtue of this fact, both K and K' are connected. You will find the following series expansion for K'(k) derived in <a href="http://www.amazon.com/AGM-Computational-Complexity-Mathematical-Monographs/dp/047131515X" rel="nofollow">Borwein's book Pi and AGM</a> Section 1.3</p> <p>$ K'(k) = \frac{2}{\pi} log \frac {4}{k} K(k) - 2 [(\frac{1}{2})^2(\frac{1}{1.2}k^2 + (\frac{1.3}{2.4})^2(\frac{1}{1.2} + \frac{1}{3.4})k^4 + (\frac{1.3.5}{2.4.6})^2(\frac{1}{1.2} + \frac{1}{3.4} + \frac{1}{5.6})k^6 $.....(infinite series) + ]</p></li> </ol> <p>You may also find Chapter 5 of <a href="http://www.amazon.com/Elliptic-Functions-Mathematical-Society-Student/dp/0521780780/ref=sr_1_1?s=books&amp;ie=UTF8&amp;qid=1280880281&amp;sr=1-1" rel="nofollow">Armitage and Eberlein's book on Elliptic Functions</a> useful.</p> <p>EDIT1: I put in the complete series expansion for K'(k). </p> <hr> <p>Regarding the computation of E(k), E(k) and K(k) are connected by the differential equation $ \frac{dK}{dk} = \frac{E - (1-k^2)K}{k(1-k^2)} $ which is how the <a href="http://mathworld.wolfram.com/LegendreRelation.html" rel="nofollow">Legendre relation</a> you mention above comes about. </p> <p>Again Borwein has the solution for this problem(buy the book!). Exercise 3 in Sec 1.4 has the formula based on the quartic AGM iteration $ E(k) = K(k)[1 - \sum_{n=0}^{\infty} 4^n [\alpha_n^4 - (\frac{\alpha_n^2+\beta_n^2}{2})^2 ] $ where</p> <ul> <li>$\alpha_n = (a_{2n})^{\frac{1}{2}} and \beta_n = (b_{2n})^{\frac{1}{2}}$</li> </ul> <p>and $ a_n, b_n$ and $c_n $ satisfy the AGM relation</p> http://mathoverflow.net/questions/34363/simultaneously-computing-a-complete-elliptic-integral-and-its-complement/34463#34463 Answer by kramyp for Simultaneously computing a complete elliptic integral and its complement kramyp 2010-08-04T01:34:25Z 2010-08-04T01:40:50Z <p>$K(b)=\int_0^{\pi/2}\frac{\mathrm{d}t}{\sqrt{\cos^2t + b\sin^2t}}$</p> <p>for K(m) -> b=1-m</p> <p>for K(1-m) -> b=m</p> http://mathoverflow.net/questions/34363/simultaneously-computing-a-complete-elliptic-integral-and-its-complement/34597#34597 Answer by J. M. for Simultaneously computing a complete elliptic integral and its complement J. M. 2010-08-05T09:21:44Z 2010-08-05T09:21:44Z <p>Short answer: yes you can. After accepting and upvoting Sandeep's question, however, the algorithm I settled on was <strong>not</strong> the algorithm he proposed.</p> <p>I managed (through a friend's assistance) to see the first two chapters of Borwein and Borwein's "Pi and the AGM". The identity that I needed, it turns out, was exactly the algorithm for the nome <code>$q$</code> that I had mentioned in my comment to his answer (i.e., a hijacked AGM iteration that computes both <code>$q=\exp(-\pi K(1-m)/K(m))$</code> and <code>$K(m)$</code> simultaneously). It is equation 2.5.14 in "Pi and the AGM":</p> <p><code>$q=\frac{m}{16}\prod_{j=0}^\infty{\left(\frac{a_j}{a_{j+1}}\right)^{2^{1-j}}}$</code></p> <p>where <code>$a_j$</code> is an iterate from the AGM iteration.</p> <p>In pseudocode, here is the algorithm I now have for computing <code>$K(m)$</code>, <code>$E(m)$</code>, and <code>$t=\ln(q(m))$</code> simultaneously:</p> <pre><code>a←(1+√(1-m))/2 c←m/(4a) t←ln(c/(4a)) s←a^2 f←1 repeat v←(a+√((a-c)(a+c)))/2 t←t+ln(a/v)/f a←v c←(c^2)/(4a) f←2f s←s-fc^2 until |c| is small enough K←π/(2a) E←KS </code></pre> <p>The code fails for <code>$m=0$</code> and <code>$m=1$</code>, but these trivial cases can be handled separately. <code>$K(1-m)$</code> is then <code>-Kt/π</code>. <code>$E(1-m)$</code> is of course then computed through the Legendre relation.</p>