Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 64
Текст из файла (страница 64)
x, y, and z must be nonnegative, and at most one can be zero. TINY must be at least 5 times the machine underflow limit,BIG at most one fifth the machine overflow limit.{float alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt;if (FMIN(FMIN(x,y),z) < 0.0 || FMIN(FMIN(x+y,x+z),y+z) < TINY ||FMAX(FMAX(x,y),z) > BIG)nrerror("invalid arguments in rf");xt=x;yt=y;zt=z;do {sqrtx=sqrt(xt);sqrty=sqrt(yt);sqrtz=sqrt(zt);alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;xt=0.25*(xt+alamb);yt=0.25*(yt+alamb);zt=0.25*(zt+alamb);ave=THIRD*(xt+yt+zt);delx=(ave-xt)/ave;dely=(ave-yt)/ave;delz=(ave-zt)/ave;} while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz)) > ERRTOL);e2=delx*dely-delz*delz;e3=delx*dely*delz;return (1.0+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave);}A value of 0.08 for the error tolerance parameter is adequate for single precision (7significant digits).
Since the error scales as 6n , we see that 0.0025 will yield double precision(16 significant digits) and require at most two or three more iterations. Since the coefficientsof the sixth-order truncation error are different for the other elliptic functions, these values forthe error tolerance should be changed to 0.04 and 0.0012 in the algorithm for RC , and 0.05 and0.0015 for RJ and RD . As well as being an algorithm in its own right for certain combinationsof elementary functions, the algorithm for RC is used repeatedly in the computation of RJ .The C implementations test the input arguments against two machine-dependent constants, TINY and BIG, to ensure that there will be no underflow or overflow during thecomputation.
We have chosen conservative values, corresponding to a machine minimumof 3 × 10−39 and a machine maximum of 1.7 × 1038 . You can always extend the range ofadmissible argument values by using the homogeneity relations (6.11.22), below.#include <math.h>#include "nrutil.h"#define ERRTOL 0.056.11 Elliptic Integrals and Jacobian Elliptic Functions#define#define#define#define#define#define#define#define265TINY 1.0e-25BIG 4.5e21C1 (3.0/14.0)C2 (1.0/6.0)C3 (9.0/22.0)C4 (3.0/26.0)C5 (0.25*C3)C6 (1.5*C4)float rd(float x, float y, float z)Computes Carlson’s elliptic integral of the second kind, RD (x, y, z). x and y must be nonnegative, and at most one can be zero.
z must be positive. TINY must be at least twice thenegative 2/3 power of the machine overflow limit. BIG must be at most 0.1 × ERRTOL timesthe negative 2/3 power of the machine underflow limit.{float alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty,sqrtz,sum,xt,yt,zt;if (FMIN(x,y) < 0.0 || FMIN(x+y,z) < TINY || FMAX(FMAX(x,y),z) > BIG)nrerror("invalid arguments in rd");xt=x;yt=y;zt=z;sum=0.0;fac=1.0;do {sqrtx=sqrt(xt);sqrty=sqrt(yt);sqrtz=sqrt(zt);alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;sum += fac/(sqrtz*(zt+alamb));fac=0.25*fac;xt=0.25*(xt+alamb);yt=0.25*(yt+alamb);zt=0.25*(zt+alamb);ave=0.2*(xt+yt+3.0*zt);delx=(ave-xt)/ave;dely=(ave-yt)/ave;delz=(ave-zt)/ave;} while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz)) > ERRTOL);ea=delx*dely;eb=delz*delz;ec=ea-eb;ed=ea-6.0*eb;ee=ed+ec+ec;return 3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*delz*ee)+delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*sqrt(ave));}#include <math.h>#include "nrutil.h"#define ERRTOL 0.05#define TINY 2.5e-13#define BIG 9.0e11#define C1 (3.0/14.0)#define C2 (1.0/3.0)#define C3 (3.0/22.0)#define C4 (3.0/26.0)#define C5 (0.75*C3)#define C6 (1.5*C4)#define C7 (0.5*C2)#define C8 (C3+C3)266Chapter 6.Special Functionsfloat rj(float x, float y, float z, float p)Computes Carlson’s elliptic integral of the third kind, RJ (x, y, z, p).
x, y, and z must benonnegative, and at most one can be zero. p must be nonzero. If p < 0, the Cauchy principalvalue is returned. TINY must be at least twice the cube root of the machine underflow limit,BIG at most one fifth the cube root of the machine overflow limit.{float rc(float x, float y);float rf(float x, float y, float z);float a,alamb,alpha,ans,ave,b,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee,fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt;if (FMIN(FMIN(x,y),z) < 0.0 || FMIN(FMIN(x+y,x+z),FMIN(y+z,fabs(p))) < TINY|| FMAX(FMAX(x,y),FMAX(z,fabs(p))) > BIG)nrerror("invalid arguments in rj");sum=0.0;fac=1.0;if (p > 0.0) {xt=x;yt=y;zt=z;pt=p;} else {xt=FMIN(FMIN(x,y),z);zt=FMAX(FMAX(x,y),z);yt=x+y+z-xt-zt;a=1.0/(yt-p);b=a*(zt-yt)*(yt-xt);pt=yt+b;rho=xt*zt/yt;tau=p*pt/yt;rcx=rc(rho,tau);}do {sqrtx=sqrt(xt);sqrty=sqrt(yt);sqrtz=sqrt(zt);alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;alpha=SQR(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz);beta=pt*SQR(pt+alamb);sum += fac*rc(alpha,beta);fac=0.25*fac;xt=0.25*(xt+alamb);yt=0.25*(yt+alamb);zt=0.25*(zt+alamb);pt=0.25*(pt+alamb);ave=0.2*(xt+yt+zt+pt+pt);delx=(ave-xt)/ave;dely=(ave-yt)/ave;delz=(ave-zt)/ave;delp=(ave-pt)/ave;} while (FMAX(FMAX(fabs(delx),fabs(dely)),FMAX(fabs(delz),fabs(delp))) > ERRTOL);ea=delx*(dely+delz)+dely*delz;eb=delx*dely*delz;ec=delp*delp;ed=ea-3.0*ec;ee=eb+2.0*delp*(ea-ec);ans=3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4))+delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave));if (p <= 0.0) ans=a*(b*ans+3.0*(rcx-rf(xt,yt,zt)));return ans;}6.11 Elliptic Integrals and Jacobian Elliptic Functions267#include <math.h>#include "nrutil.h"#define ERRTOL 0.04#define TINY 1.69e-38#define SQRTNY 1.3e-19#define BIG 3.e37#define TNBG (TINY*BIG)#define COMP1 (2.236/SQRTNY)#define COMP2 (TNBG*TNBG/25.0)#define THIRD (1.0/3.0)#define C1 0.3#define C2 (1.0/7.0)#define C3 0.375#define C4 (9.0/22.0)float rc(float x, float y)Computes Carlson’s degenerate elliptic integral, RC (x, y).
x must be nonnegative and y mustbe nonzero. If y < 0, the Cauchy principal value is returned. TINY must be at least 5 timesthe machine underflow limit, BIG at most one fifth the machine maximum overflow limit.{float alamb,ave,s,w,xt,yt;if (x < 0.0 || y == 0.0 || (x+fabs(y)) < TINY || (x+fabs(y)) > BIG ||(y<-COMP1 && x > 0.0 && x < COMP2))nrerror("invalid arguments in rc");if (y > 0.0) {xt=x;yt=y;w=1.0;} else {xt=x-y;yt = -y;w=sqrt(x)/sqrt(xt);}do {alamb=2.0*sqrt(xt)*sqrt(yt)+yt;xt=0.25*(xt+alamb);yt=0.25*(yt+alamb);ave=THIRD*(xt+yt+yt);s=(yt-ave)/ave;} while (fabs(s) > ERRTOL);return w*(1.0+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave);}At times you may want to express your answer in Legendre’s notation. Alternatively, you may be given results in that notation and need to compute their valueswith the programs given above.
It is a simple matter to transform back and forth.The Legendre elliptic integral of the 1st kind is defined as&F (φ, k) ≡0φ(dθ1 − k 2 sin2 θ(6.11.17)The complete elliptic integral of the 1st kind is given byK(k) ≡ F (π/2, k)(6.11.18)In terms of RF ,F (φ, k) = sin φRF (cos2 φ, 1 − k 2 sin2 φ, 1)K(k) = RF (0, 1 − k 2 , 1)(6.11.19)268Chapter 6.Special FunctionsThe Legendre elliptic integral of the 2nd kind and the complete elliptic integral ofthe 2nd kind are given by&φE(φ, k) ≡(1 − k 2 sin2 θ dθ0= sin φRF (cos2 φ, 1 − k 2 sin2 φ, 1)(6.11.20)− 13 k 2 sin3 φRD (cos2 φ, 1 − k 2 sin2 φ, 1)E(k) ≡ E(π/2, k) = RF (0, 1 − k 2 , 1) − 13 k 2 RD (0, 1 − k 2 , 1)Finally, the Legendre elliptic integral of the 3rd kind is&φΠ(φ, n, k) ≡0dθ((1 + n sin θ) 1 − k 2 sin2 θ2= sin φRF (cos2 φ, 1 − k 2 sin2 φ, 1)(6.11.21)− 13 n sin3 φRJ (cos2 φ, 1 − k 2 sin2 φ, 1, 1 + n sin2 φ)(Note that this sign convention for n is opposite that of Abramowitz and Stegun [12],and that their sin α is our k.)#include <math.h>#include "nrutil.h"float ellf(float phi, float ak)Legendre elliptic integral of the 1st kind F (φ, k), evaluated using Carlson’s function RF .
Theargument ranges are 0 ≤ φ ≤ π/2, 0 ≤ k sin φ ≤ 1.{float rf(float x, float y, float z);float s;s=sin(phi);return s*rf(SQR(cos(phi)),(1.0-s*ak)*(1.0+s*ak),1.0);}#include <math.h>#include "nrutil.h"float elle(float phi, float ak)Legendre elliptic integral of the 2nd kind E(φ, k), evaluated using Carlson’s functions RD andRF . The argument ranges are 0 ≤ φ ≤ π/2, 0 ≤ k sin φ ≤ 1.{float rd(float x, float y, float z);float rf(float x, float y, float z);float cc,q,s;s=sin(phi);cc=SQR(cos(phi));q=(1.0-s*ak)*(1.0+s*ak);return s*(rf(cc,q,1.0)-(SQR(s*ak))*rd(cc,q,1.0)/3.0);}6.11 Elliptic Integrals and Jacobian Elliptic Functions269#include <math.h>#include "nrutil.h"float ellpi(float phi, float en, float ak)Legendre elliptic integral of the 3rd kind Π(φ, n, k), evaluated using Carlson’s functions RJ andRF . (Note that the sign convention on n is opposite that of Abramowitz and Stegun.) Theranges of φ and k are 0 ≤ φ ≤ π/2, 0 ≤ k sin φ ≤ 1.{float rf(float x, float y, float z);float rj(float x, float y, float z, float p);float cc,enss,q,s;s=sin(phi);enss=en*s*s;cc=SQR(cos(phi));q=(1.0-s*ak)*(1.0+s*ak);return s*(rf(cc,q,1.0)-enss*rj(cc,q,1.0,1.0+enss)/3.0);}Carlson’s functions are homogeneous of degree − 12 and − 32 , soRF (λx, λy, λz) = λ−1/2 RF (x, y, z)RJ (λx, λy, λz, λp) = λ−3/2 RJ (x, y, z, p)(6.11.22)Thus to express a Carlson function in Legendre’s notation, permute the first threearguments into ascending order, use homogeneity to scale the third argument to be1, and then use equations (6.11.19)–(6.11.21).Jacobian Elliptic FunctionsThe Jacobian elliptic function sn is defined as follows: instead of consideringthe elliptic integralu(y, k) ≡ u = F (φ, k)(6.11.23)consider the inverse functiony = sin φ = sn(u, k)(6.11.24)Equivalently,& snu=0dy(2(1 − y )(1 − k 2 y2 )(6.11.25)When k = 0, sn is just sin.
The functions cn and dn are defined by the relationssn2 + cn2 = 1,k 2 sn2 + dn2 = 1(6.11.26)The routine given below actually takes mc ≡ kc2 = 1 − k 2 as an input parameter.It also computes all three functions sn, cn, and dn since computing all three is noharder than computing any one of them. For a description of the method, see [8].270Chapter 6.#include <math.h>#define CA 0.0003Special FunctionsThe accuracy is the square of CA.void sncndn(float uu, float emmc, float *sn, float *cn, float *dn)Returns the Jacobian elliptic functions sn(u, kc ), cn(u, kc ), and dn(u, kc ). Here uu = u, whileemmc = kc2 .{float a,b,c,d,emc,u;float em[14],en[14];int i,ii,l,bo;emc=emmc;u=uu;if (emc) {bo=(emc < 0.0);if (bo) {d=1.0-emc;emc /= -1.0/d;u *= (d=sqrt(d));}a=1.0;*dn=1.0;for (i=1;i<=13;i++) {l=i;em[i]=a;en[i]=(emc=sqrt(emc));c=0.5*(a+emc);if (fabs(a-emc) <= CA*a) break;emc *= a;a=c;}u *= c;*sn=sin(u);*cn=cos(u);if (*sn) {a=(*cn)/(*sn);c *= a;for (ii=l;ii>=1;ii--) {b=em[ii];a *= c;c *= (*dn);*dn=(en[ii]+a)/(b+a);a=c/b;}a=1.0/sqrt(c*c+1.0);*sn=(*sn >= 0.0 ? a : -a);*cn=c*(*sn);}if (bo) {a=(*dn);*dn=(*cn);*cn=a;*sn /= d;}} else {*cn=1.0/cosh(u);*dn=(*cn);*sn=tanh(u);}}6.12 Hypergeometric Functions271CITED REFERENCES AND FURTHER READING:Erdélyi, A., Magnus, W., Oberhettinger, F., and Tricomi, F.G.
1953, Higher TranscendentalFunctions, Vol. II, (New York: McGraw-Hill). [1]Gradshteyn, I.S., and Ryzhik, I.W. 1980, Table of Integrals, Series, and Products (New York:Academic Press). [2]Carlson, B.C. 1977, SIAM Journal on Mathematical Analysis, vol. 8, pp. 231–242. [3]Carlson, B.C. 1987, Mathematics of Computation, vol.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
