forum.alglib.net
http://forum.alglib.net/

Error in elliptical integral formula
http://forum.alglib.net/viewtopic.php?f=2&t=566		 Page 1 of 1
Author:  mikefee [ Thu Apr 19, 2012 2:44 am ]
Post subject:  Error in elliptical integral formula
I just tried using the Delphi version of the following functions (from the unit elliptic.pas)

function EllipticIntegralK(m : Double):Double;
function EllipticIntegralE(m : Double):Double;

to calculate the complete elliptic integrals of the first and second type respectively. It appears that the functions do not return the correct numerical values. I used the following procedure to compare Alglib output with an alternate iterative solution...

procedure TestElliptic;
var
i : integer;
k, AlglibK, AlglibE, AltK, AltE : real;
begin
Form1.memo1.lines.clear;
Form1.memo1.lines.add('---k-------AlglibK---------AlglibE-----------AltK-----------AltE');
for i := 1 to 19 do
begin
k := i/20.0;
AlglibK := EllipticIntegralK(k);
AlglibE := EllipticIntegralE(k);
AltK := AltEllipticIntegralK(k);
AltE := AltEllipticIntegralE(k);
Form1.memo1.lines.add(Format('%5.2f : %12.9f %12.9f : %12.9f %12.9f',
[k,AlglibK,AlglibE,AltK,AltE]));
end;
end;

...I get the following results:

---k-------AlglibK---------AlglibE-----------AltK-----------AltE
0.05 : 1.591003454 1.550973352 : 1.571779457 1.569814118
0.10 : 1.612441349 1.530757637 : 1.574745562 1.566861942
0.15 : 1.635256732 1.510121832 : 1.579745661 1.561922968
0.20 : 1.659623599 1.489035058 : 1.586867847 1.554968546
0.25 : 1.685750355 1.467462209 : 1.596242222 1.545957256
0.30 : 1.713889448 1.445363064 : 1.608048620 1.534833465
0.35 : 1.744350597 1.422691133 : 1.622528103 1.521525269
0.40 : 1.777519371 1.399392139 : 1.639999866 1.505941612
0.45 : 1.813883937 1.375401972 : 1.660886238 1.487968272
0.50 : 1.854074677 1.350643881 : 1.685750355 1.467462209
0.55 : 1.898924910 1.325024498 : 1.715354496 1.444243488
0.60 : 1.949567750 1.298428035 : 1.750753803 1.418083394
0.65 : 2.007598398 1.270707480 : 1.793454100 1.388686356
0.70 : 2.075363135 1.241670568 : 1.845693998 1.355661136
0.75 : 2.156515647 1.211056028 : 1.910989781 1.318472108
0.80 : 2.257205327 1.178489924 : 1.995302778 1.276349943
0.85 : 2.389016486 1.143395792 : 2.109935468 1.228108314
0.90 : 2.578092113 1.104774733 : 2.280549138 1.171697053
0.95 : 2.908337248 1.060473728 : 2.590011231 1.102721648

Comparison of Alglib output with an alternate algorithm plus published tables and other resources (for example the calculator at http://keisan.casio.com/has10/SpecExec.cgi) suggests that the Alglib results (on the left of the table) are incorrect.

Mike
Page 1 of 1 All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/