I have noticed a rather significant discrepancy between the values of the elliptic integrals returned by ellipke() and those published in a number of texts. I am not sure if this is me not understanding the inputs or if it is genuinely wrong. I have checked these values against two texts as well as my own calculation based on Simpson's rule.
If anyone could offer some advice or tell me what I am doing wrong, it would be VERY appreciated.
_________________________________________________________________
Reference 1: Lewis,'Vortex Element Method for Fluid Dynamic Analysis of Engineering Systems' 1991 pg 519-520.
Reference 2: Liepman, 'Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. pg 610-611. <http://www.convertit.com/Go/GovCon/Reference/AMS55.ASP?Res=150&Page=610&Submit=Go> _________________________________________________________________
For alpha=45 deg, k=sin(alpha)=0.707107
Reference 1 Gives: K=1.854074677301372 and E=1.350643881047676
Reference 2 Gives: K=1.854075 and E=1.350644
[K,E]=ellipke(0.707107); Gives: K=2.085974029127680 and E=1.237422393581531
Integrating Via Simpson's Rule Gives: K=1.854075629303571 and E=1.350644357747139
_________________________________________________________________
My Simpsons Rule Code:
function [ K,E ] = FunEllipticIntegrals( phi )%EllipticalIntegral1stKind Calculates the value of the elliptic integral of
%the first kind using simpsons rule.
tol=10^-6;k=sind(phi);i=1;for a=0:tol:pi/2 b=a+tol; [FunK1] = FunFirstkind(k,a); [FunK2] = FunFirstkind(k,(a+b)/2); [FunK3] = FunFirstkind(k,b); FunK(i) = (b-a)/6*(FunK1+4*FunK2+FunK3); [FunE1] = FunSecondkind(k,a); [FunE2] = FunSecondkind(k,(a+b)/2); [FunE3] = FunSecondkind(k,b); FunE(i) = (b-a)/6*(FunE1+4*FunE2+FunE3); i=i+1;endK=sum(FunK);E=sum(FunE);endfunction [ FunK ] = FunFirstkind( k,alpha )%FunFirstkind Elliptic Integral of 1st Kind Fuction
FunK=1/sqrt(1-k^2*sin(alpha)^2);endfunction [ FunE ] = FunSecondkind( k,alpha )%FunFirstkind Elliptic Integral of 2nd Kind Fuction
FunE=sqrt(1-k^2*sin(alpha)^2);end
Best Answer