[Math] Heuman Lambda Function in MATLAB

elliptic integralsintegrationMATLABspecial functions

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!

Best Answer

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

Related Question