home *** CD-ROM | disk | FTP | other *** search
/ Chip 2002 June / Chip_2002-06_cd1.bin / zkuste / derive / download / Setup.exe / %MAINDIR% / Users / Carlson.mth < prev    next >
Encoding:
Text File  |  2002-05-15  |  6.5 KB  |  164 lines
  1. "Carlson.mth Carlson elliptic integrals."
  2.  
  3. "These numerical algroithms are based on Numerical Algorithms 10(1995)13-26."
  4.  
  5. "B.C. Carlson Numerical computation of real or complex elliptic integrals"
  6.  
  7. "They have been updated to the algorithms in Journal of Computational"
  8.  
  9. "and Applied Mathematics 118 (2000) 71-85"
  10.  
  11. "Reduction theorms for elliptic integrands with the square root of"
  12.  
  13. "two quadritic factors"
  14.  
  15. "B.C. Carlson, James FitzSimons"
  16.  
  17. "Algorithm for R_F"
  18.  
  19. "page 16 2.3"
  20.  
  21. "page 16 2.4"
  22.  
  23. "Iterative definition of R_F"
  24.  
  25. R_FS9(e2,e3):=195*e2^8/32768-429*e2^7/59392+e2^6*(3003*e3/63488+231/25600)-e2~
  26. ^5*(819*e3^2+308*e3+72)/6144+e2^4*(3465*e3^2/29696+315*e3/5888+35/2176)-e2^3*~
  27. (1426425*e3^3+964782*e3^2+564200*e3+235600)/9800960+e2^2*(105*e3^4/1024+35*e3~
  28. ^3/384+5*e3^2/64+e3/16+1/24)-e2*(6774075*e3^4+7592200*e3^3+8804400*e3^2+10885~
  29. 440*e3+15965312)/159653120+63*e3^5/7936+7*e3^4/640+5*e3^3/304+3*e3^2/104+e3/1~
  30. 4
  31.  
  32. R_F(xx,yy,zz,v,v1,q,a,n,lambda,x,y,z,e2,e3):=PROG(v:=[xx,yy,zz],q:=VECTOR(ABS~
  33. (x),x,v),a:=MIN(q),q:=MAX(q),n:=CEILING((LN(IF(a>0,q/a,q))+LN(10)*PrecisionDi~
  34. gits/16)/LN(4))+4,LOOP(v1:=VECTOR(SQRT(i_),i_,v),lambda:=v1 SUB 1*v1 SUB 2+v1~
  35.  SUB 2*v1 SUB 3+v1 SUB 3*v1 SUB 1,v:=VECTOR(i_+lambda,i_,v)/4,n:-1,IF(n=0,exi~
  36. t)),a:=AVERAGE(v),x:=1-v SUB 1/a,y:=1-v SUB 2/a,z:=-x-y,e2:=x*y-z^2,e3:=x*y*z~
  37. ,a:=1/SQRT(a),a*R_FS9(e2,e3)+a)
  38.  
  39. "Algorithm for R_C"
  40.  
  41. "page 17 2.10"
  42.  
  43. "Iterative definition of R_C"
  44.  
  45. "Add a special case when y is real and negative."
  46.  
  47. "page 18 2.14"
  48.  
  49. R_C_AUX2(x):=SUM((-x)^i_/(2*i_+1),i_,1,1+ABS(PrecisionDigits*LN(10)/(-LN(ABS(~
  50. x)+0.01))))
  51.  
  52. R_C_AUX5(xx,yy,x,y,lambda):=PROG(x:=xx,y:=yy,LOOP(IF(ABS((y-x)/x)<1/64,RETURN~
  53. ((R_C_AUX2((y-x)/x)+1)/SQRT(x))),lambda:=2*SQRT(x)*SQRT(y)+y,x:=(x+lambda)/4,~
  54. y:=(y+lambda)/4))
  55.  
  56. R_C(x,y):=IF(IM(y)=0 AND RE(y)<0,IF(x=0,0,SQRT(x/(x-y))*R_C_AUX5(x-y,-y)),IF(~
  57. x=0,pi/(2*SQRT(y)),R_C_AUX5(x,y)))
  58.  
  59. "Algorithm for R_J"
  60.  
  61. "page 18 2.17"
  62.  
  63. "page 19 2.21"
  64.  
  65. "Iterative definition of R_J"
  66.  
  67. S_M2(x):=IF(ABS(x)<1/64,SUM((-x/4)^i_*COMB(2*i_,i_)/(1-2*i_),i_,1,1+ABS(Preci~
  68. sionDigits*LN(10)/(-LN(ABS(x)+0.01)))),SQRT(1+x)-1)
  69.  
  70. R_JS9(e2,e3,e4,e5):=3861*e2^8/229376-1287*e2^7/63488+e2^6*(273*e3/2048-1287*e~
  71. 4/10240+77/3072)-9*e2^5*(8870433*e3^2+3294060*e3-4060*(759*e4-713*e5-186))/21~
  72. 1732480+e2^4*(10395*e3^2/31744+9*e3*(3*(55*e5+14)-175*e4)/2560+297*e4^2/1024-~
  73. 35*e4/256+945*e5/7424+105/2432)-e2^3*(119409675*e3^3+1137235*e3^2*(70-297*e4)~
  74. -25300*e3*(5859*e4-5481*e5-1798)+58*(1195425*e4^2-32550*e4*(69*e5+22)+7843*(1~
  75. 35*e5^2+84*e5+40)))/291132160+e2^2*(297*e3^4/1024+945*e3^3/3712+315*e3^2*(31*~
  76. (69*e5+22)-2277*e4)/1003904+e3*(945*e4^2/1408-9*e4*(45*e5+14)/320+35*e5/96+45~
  77. /272)-27*e4^3/128+35*e4^2/192-45*e4*(133*e5+58)/17632+315*e5^2/1984+15*e5/112~
  78. +9/88)-e2*(3056930051175*e3^4-21879326196*e3^3*(525*e4-495*e5-154)+121331210*~
  79. e3^2*(133893*e4^2-77140*e4+540*(133*e5+58))+116280*e3*(74939865*e4^2-953810*e~
  80. 4*(147*e5+62)+6293*(10465*e5^2+8580*e5+6072))-57304*(47418525*e4^3-830025*e4^~
  81. 2*(161*e5+66)+8399853*e4*(15*e5^2+12*e5+8)-46665850*e5^2-59293080*e5-95998320~
  82. ))/25671742736640+63*e3^5/2816+e3^4*(70-243*e4)/2304+5*e3^3*(29*(147*e5+62)-4~
  83. 557*e4)/201376+e3^2*(315*e4^2/1984-15*e4*(161*e5+66)/8096+9*e5^2/64+9*e5/80+3~
  84. /40)-e3*(4917675*e4^3-2781999*e4^2*(5*e5+2)+121220*e4*(85*e5+54)-1870*(2565*e~
  85. 5^2+3132*e5+4408))/49457760+3*e4^4/128-5*e4^3/144+e4^2*(45*e5/464+9/152)-3*e4~
  86. *(1155*e5^2+1364*e5+1736)/38192+5*e5^3/176+9*e5^2/184+3*e5/26
  87.  
  88. R_J_AUX6(v,q,a,n,v1,x,y,z,p,lambda,alpha,sm,delta,m,a1,e2,e3,e4,e5):=PROG(q:=~
  89. VECTOR(ABS(x),x,v),a:=MIN(q),q:=MAX(q),n:=CEILING((LN(IF(a>0,q/a,q))+LN(10)*P~
  90. recisionDigits/16)/LN(2))+4,x:=SQRT(v SUB 1),y:=SQRT(v SUB 2),z:=SQRT(v SUB 3~
  91. ),p:=v SUB 4,lambda:=x*y+y*z+z*x,v1:=VECTOR(i_+lambda,i_,v)/4,alpha:=p*(x+y+z~
  92. )+x*y*z,delta:=PRODUCT(p-v SUB i_,i_,1,3),sm:=alpha+alpha*S_M2(delta/alpha^2)~
  93. /2,delta:=delta/4,m:=1/4,LOOP(x:=SQRT(v1 SUB 1),y:=SQRT(v1 SUB 2),z:=SQRT(v1 ~
  94. SUB 3),p:=SQRT(v1 SUB 4),lambda:=x*y+y*z+z*x,dm:=(p+x)*(p+y)*(p+z),rm:=2*sm+s~
  95. m*S_M2(delta/sm^2),sm:=(dm*rm-m*delta)/(2*(dm+m*rm)),delta:=delta/4,m:=m/4,v1~
  96. :=VECTOR(i_+lambda,i_,v1)/4,n:-1,IF(n<=0,exit)),x:=v1 SUB 1,y:=v1 SUB 2,z:=v1~
  97.  SUB 3,p:=v1 SUB 4,a:=(x+y+z+2*p)/5,x:=1-x/a,y:=1-y/a,z:=1-z/a,p:=(-x-y-z)/2,~
  98. e2:=-3*p^2+z*(y+x)+x*y,e3:=2*p*e2+4*p^3+x*y*z,e4:=p*(p*e2+3*p^3+2*x*y*z),e5:=~
  99. x*y*z*p^2,a:=1/a^(3/2),y:=delta/sm^2,m*(a*R_JS9(e2,e3,e4,e5)+a)+3/sm*IF(ABS(y~
  100. )<1/64,R_C_AUX2(y)+1,R_C_AUX5(1,1+y)))
  101.  
  102. "Add a special case when p is real and negative."
  103.  
  104. "page 19 2.26"
  105.  
  106. SORT_RJ(v):=IF(RE(v SUB 1)>RE(v SUB 2),SORT_RJ([v SUB 2,v SUB 1,v SUB 3,v SUB~
  107.  4]),IF(RE(v SUB 2)>RE(v SUB 3),SORT_RJ([v SUB 1,v SUB 3,v SUB 2,v SUB 4]),[v~
  108.  SUB 1,v SUB 2,v SUB 3,-v SUB 4,(v SUB 1*(v SUB 2-v SUB 3)+v SUB 2*(v SUB 3-v~
  109.  SUB 4))/(v SUB 2-v SUB 4)]))
  110.  
  111. R_J(x,y,z,p,v):=IF(IM(p)=0 AND RE(p)<0,PROG(v:=SORT_RJ([x,y,z,p]),((v SUB 5-v~
  112.  SUB 2)*R_J_AUX6([v SUB 1,v SUB 2,v SUB 3,v SUB 5])-3*R_F(v SUB 1,v SUB 2,v S~
  113. UB 3)+3*SQRT(v SUB 1*v SUB 2*v SUB 3/(v SUB 1*v SUB 3+v SUB 5*v SUB 4))*R_C(v~
  114.  SUB 1*v SUB 3+v SUB 5*v SUB 4,v SUB 5*v SUB 4))/(v SUB 2+v SUB 4)),R_J_AUX6(~
  115. [x,y,z,p]))
  116.  
  117. "Algorithm for R_D"
  118.  
  119. "page 20 2.28"
  120.  
  121. "page 20 2.27"
  122.  
  123. "Iterative definition of R_D"
  124.  
  125. R_D(xx,yy,zz,v,v1,v2,q,a,x,y,z,lambda,s,f,e2,e3,e4,e5):=PROG(v:=[xx,yy,zz],q:~
  126. =VECTOR(ABS(x),x,v),a:=MIN(q),q:=MAX(q),n:=CEILING((LN(IF(a>0,q/a,q))+LN(10)*~
  127. PrecisionDigits/16)/LN(4))+4,s:=0,f:=1,v1:=v,LOOP(x:=SQRT(v1 SUB 1),y:=SQRT(v~
  128. 1 SUB 2),z:=SQRT(v1 SUB 3),lambda:=x*y+y*z+z*x,s:+f*3/(z*(v1 SUB 3+lambda)),v~
  129. 1:=VECTOR(i_+lambda,i_,v1)/4,f:/4,n:-1,IF(n=0,exit)),x:=v1 SUB 1,y:=v1 SUB 2,~
  130. z:=v1 SUB 3,a:=(x+y+3*z)/5,x:=1-x/a,y:=1-y/a,z:=(-x-y)/3,e2:=x*y-6*z^2,e3:=(3~
  131. *x*y-8*z^2)*z,e4:=3*(x*y-z^2)*z^2,e5:=x*y*z^3,a:=1/a^(3/2),f*(a*R_JS9(e2,e3,e~
  132. 4,e5)+a)+s)
  133.  
  134. "Algorithm for R_G"
  135.  
  136. "page 15 1.7"
  137.  
  138. R_G(x,y,z):=(z*R_F(x,y,z)-(x-z)*(y-z)*R_D(x,y,z)/3+SQRT(x*y/z))/2
  139.  
  140. "page 24 4. Other integrals"
  141.  
  142. "These are the Legendre elliptic integrals."
  143.  
  144. LEGENDRE_F(phi,m):=R_F(CSC(phi)^2-1,CSC(phi)^2-m,CSC(phi)^2)
  145.  
  146. LEGENDRE_E(phi,m):=LEGENDRE_F(phi,m)-m/3*R_D(CSC(phi)^2-1,CSC(phi)^2-m,CSC(ph~
  147. i)^2)
  148.  
  149. LEGENDRE_PI(phi,m,n):=LEGENDRE_F(phi,m)+n/3*R_J(CSC(phi)^2-1,CSC(phi)^2-m,CSC~
  150. (phi)^2,CSC(phi)^2-n)
  151.  
  152. ELLIPTICF(x,k):=R_F(1/x^2-1,1/x^2-k^2,1/x^2)
  153.  
  154. ELLIPTICE(x,k):=ELLIPTICF(x,k)-k^2/3*R_D(1/x^2-1,1/x^2-k^2,1/x^2)
  155.  
  156. ELLIPTICPI(x,n,k):=ELLIPTICF(x,k)+n/3*R_J(1/x^2-1,1/x^2-k^2,1/x^2,1/x^2-n)
  157.  
  158. ELLIPTICK(k):=R_F(0,1-k^2,1)
  159.  
  160. ELLIPTICEK(k):=ELLIPTICK(k)-k^2/3*R_D(0,1-k^2,1)
  161.  
  162. ELLIPTICPIK(n,k):=ELLIPTICK(k)+n/3*R_J(0,1-k^2,1,1-n)
  163.  
  164.