Heuman Lambda Function in MATLAB

I'm trying to implement the Heuman Lambda function in MATLAB, but i'm having problems getting a correct answer.

The Heuman Lambda Function is:

$$\Lambda_0(\beta,k) = \frac{2}{pi}[E(k)F(\beta,k') + K(k)E(\beta,k') - K(k)F(\beta,k')]$$

Where E(k) and K(k) are the complete elliptic integrals of the first and second kind, and E($\beta$,k) and F($\beta$,k) are the incomplete elliptic integrals of the first and second kind.

K(k) = complete elliptic integral of the first kind

$$K(k) = \int_0^{2\pi} \frac {1}{\sqrt{1 - k^2 sin^2 \theta} }d\theta$$

E(k) = complete elliptic integral of the second kind

$$E(k) = \int_0^{2\pi} \sqrt{1 - k^2 sin^2 \theta }d\theta$$

F($\beta$,k) = incomplete elliptic integral of the first kind

$$F(\beta,k) = \int_0^{\beta} \frac {1}{\sqrt{1 - k^2 sin^2 \theta} }d\theta$$

E($\beta$,k) = incomplete elliptic integral of the second kind

$$E(\beta,k) = \int_0^{\beta} {\sqrt{1 - k^2 sin^2 \theta} }d\theta$$

Here is my attempt so far

function [ HL ] = Heuman_Lambda( b,k )

b = b*pi/180;    %see note 1

k = sin(k*pi/180)^2;    %see note 2

kdash = (1-(k^2))^0.5;    %see note 3

[K,E] = ellipke(k);

incF = ellipticF(b,kdash);

incE = ellipticE(b,kdash);

HL = 2/pi * (E*incF + K*incE - K*incF );


note 1: This is because MATLAB works in radians, but the table of answers I have uses degrees for $\beta$

note 2: This is because in the MATLAB implementation of the 1st and 2nd complete elliptic integrals the equation is slightly different. See MATLAB ellpipke at the bottom for an explanation. Also see a previous question I have asked on the subject where it was explained to me.

note 3 :This is because in This Document on page 256, just under equation 17, it says that this will compute k'

I have tried many variations of this code, but I cannot get an accurate answer. I am certain that the individual elliptic integral functions are correct for k and $\beta$, as I have individually tested all four of them and compared them to a table of answers I have. I do not have a table of answers for k' however so this may be where I am going wrong.

Here is the table of answers I have for the Heuman Lambda function:

(sin^-1)k                b
-     |    0    |    5    |    10
----------------------------------
0     |    0    | 0.08716 | 0.17364
30    |    0    | 0.08143 | 0.16225
60    |    0    | 0.06723 | 0.13413

(sin^-1)k and b are both in degrees


These do not agree with the answers I am getting for my function. Does anybody know where I am going wrong?

Thanks!

• Your k is what the theory calls $k^2$, so your kdash should use k instead of k^2 (the latter would represent $k^4$) – ccorn Jan 5 '14 at 21:00
• And since Matlab wants to be given $k^2$ in ellipke and probably in ellipticF and ellipticE as well, it probably also wants to be given $k'^2=1-k^2$ which is 1-k in your notation. That is, no explicit squareroot seems to be needed. – ccorn Jan 5 '14 at 21:08
• Thank you so much! The function now works perfectly. You have no idea how long i've been struggling with this. God I love this website. – Blue7 Jan 5 '14 at 22:36
• For reference. – J. M. is a poor mathematician May 19 '17 at 8:40

I'll post an answer so that I can close this question.

The code to implement the Heuman Lambda Function in MATLAB is:

function [ HL ] = Heuman_Lambda( b,k )

b = b*pi/180;

k = sin(k*pi/180)^2;

kdash = (1-k);

[K,E] = ellipke(k);

incF = ellipticF(b,kdash);

incE = ellipticE(b,kdash);

HL = 2/pi * (E*incF + K*incE - K*incF );


Be aware that the elliptic integrals ellipicF and ellipticE do not work straight away. You first have to create a matlab function for these calling them from the MUPAD. To do this:

function [ F ] = ellipticF( phi,k )

y = feval(symengine,'ellipticF',phi,k);

F = double(y);

end


And do the same for ellipticE