数值计算程序大放送-特殊函数

选自<<徐世良数值计算程序集(C)>>

所有的函数声明部分

//
//功能:计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷
//描述:gamma[x]=Integrate[exp[-t]*t^(x-1),{t,0,∞}]
//调用:
double lagam(double x);
//
//功能:不完全伽马函数
//描述:gamma[a,x]=P[a,x]/gamma[x]
//描述:P[a,x]=Integrate[exp[-t]*t^(a-1),{t,0,x}]
//参数:a-参数
//调用:lagam(x)函数
double lbgam(double a,double x);
//
//功能:误差函数
//描述:erf[x]=gamma[0.5,x^2]
//描述:erf[x]=2/sqrt[pi]*Integrate[exp[-t^2],{t,0,x}]
//调用:lagam(),lbgam()
double lcerf(double x);
//
//功能:第一类整数阶贝塞尔函数
//参数:n-阶数
//调用:
double ldbesl(int n,double x);
//
//功能:第二类整数阶贝塞尔函数
//参数:n-阶数
//调用:ldbesl()
double lebesl(int n,double x);
//
//功能:变型第一类整数阶贝塞尔函数
//参数:n-阶数
//调用:
double lfbesl(int n,double x);
//
//功能:变型第二类整数阶贝塞尔函数
//参数:n-阶数
//调用:lfbesl();
double lgbesl(int n,double x);
//
//功能:不完全贝塔(beta)函数
//描述:Bx[a,b]=Integrate[t^(a-1)*(1-t)^(b-1),{t,0,x}]/B[a,b]
//描述:B[a,b]=gamma[a]*gamma[b]/gamma[a+b]
//参数:a-参数,b-参数
//调用:lagam();
double lhbeta(double a,double b,double x);
//
//功能:正态分布函数
//参数:a-均值,b-方差
//调用:lcerf(),lagam(),lbgam();
double ligas(double a,double d,double x);
//
//功能:t-分布函数
//参数:n-自由度
//调用:lhbeta(),lagam();
double ljstd(double t,int n);
//
//功能:X^2-分布函数
//参数:n-自由度
//调用:lbgam(),lagam();
double lkchi(double x,int n);
//
//功能:F-分布函数
//参数:n1-自由度,n2-自由度
//调用:lhbeta(),lagam();
double llf(double f,int n1,int n2);
//
//功能:正弦积分
//参数:
//调用:
double lmsi(double x);
//
//功能:余弦积分
//参数:
//调用:
double lnci(double x);
//
//功能:指数积分
//参数:
//调用:
double loei(double x);
//
//功能:第一类椭圆积分
//参数:
//调用:
double lpfk(double k,double f);
//
//功能:第二类椭圆积分
//参数:
//调用:
double lqek(double k,double f);

 

#include "stdio.h"
#include "math.h"
//功能:计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷
//描述:Integrate[exp[-t]*t^(x-1),{t,0,∞}]
//调用:
double lagam(double x)
{
int i;
double y,t,s,u;
static double a[11]={ 0.0000677106,-0.0003442342,
  0.0015397681,-0.0024467480,0.0109736958,
  -0.0002109075,0.0742379071,0.0815782188,
  0.4118402518,0.4227843370,1.0};
if (x<=0.0)
{
  printf("err**x<=0!/n");
  return(-1.0);
}
y=x;
if (y<=1.0)
{
  t=1.0/(y*(y+1.0));
  y=y+2.0;
}
else if (y<=2.0)
{
  t=1.0/y;
  y=y+1.0;
}
else if (y<=3.0)
{
  t=1.0;
}
else
{
  t=1.0;
  while (y>3.0)
  {
   y=y-1.0;
   t=t*y;
  }
}
s=a[0];
u=y-2.0;
for (i=1; i<=10; i++)
{
  s=s*u+a[i];
}
s=s*t;
return(s);
} //功能:不完全伽马函数
//描述:gamma[a,x]=P[a,x]/gamma[x]
//描述:P[a,x]=Integrate[exp[-t]*t^(a-1),{t,0,x}]
//参数:a-参数
//调用:lagam(x)函数
double lbgam(double a,double x)
double a,x;
{
int n;
    double p,q,d,s,s1,p0,q0,p1,q1,qq;
    if ((a<=0.0)||(x<0.0))
{
  if (a<=0.0)
  {
   printf("err**a<=0!/n");
  }
        if (x<0.0)
  {
   printf("err**x<0!/n");
  }
        return(-1.0);
}
    if (x+1.0==1.0)
{
  return(0.0);
}
    if (x>1.0e+35)
{
  return(1.0);
}
    q=log(x);
q=a*q;
qq=exp(q);
    if (x<1.0+a)
{
  p=a;
  d=1.0/a;
  s=d;
        for (n=1; n<=100; n++)
  {
   p=1.0+p;
   d=d*x/p;
   s=s+d;
   if (fabs(d)<fabs(s)*1.0e-07)
   {
    s=s*exp(-x)*qq/lagam(a);
                return(s);
   }
  }
}
    else
{
  s=1.0/x;
  p0=0.0;
  p1=1.0;
  q0=1.0;
  q1=x;
        for (n=1; n<=100; n++)
  {
   p0=p1+(n-a)*p0;
   q0=q1+(n-a)*q0;
            p=x*p0+n*p1;
   q=x*q0+n*q1;
            if (fabs(q)+1.0!=1.0)
   {
    s1=p/q;
    p1=p;
    q1=q;
                if (fabs((s1-s)/s1)<1.0e-07)
    {
     s=s1*exp(-x)*qq/lagam(a);
                    return(1.0-s);
    }
                s=s1;
   }
            p1=p;
   q1=q;
  }
}
    printf("a too large !/n");
    s=1.0-s*exp(-x)*qq/lagam(a);
    return(s);
} //功能:误差函数
//描述:erf[x]=gamma[0.5,x^2]
//描述:erf[x]=2/sqrt[pi]*Integrate[exp[-t^2],{t,0,x}]
//调用:lagam(),lbgam()
double lcerf(double x)
{
double y;
    if (x>=0.0)
{
  y=lbgam(0.5,x*x);
}
    else
{
  y=-lbgam(0.5,x*x);
}
    return(y);
} //功能:第一类整数阶贝塞尔函数
//参数:n-阶数
//调用:
double ldbesl(int n,double x)
{
int i,m;
    double t,y,z,p,q,s,b0,b1;
    static double a[6]={ 57568490574.0,-13362590354.0,
  651619640.7,-11214424.18,77392.33017,-184.9052456};
    static double b[6]={ 57568490411.0,1029532985.0,
  9494680.718,59272.64853,267.8532712,1.0};
    static double c[6]={ 72362614232.0,-7895059235.0,
  242396853.1,-2972611.439,15704.4826,-30.16036606};
    static double d[6]={ 144725228443.0,2300535178.0,
  18583304.74,99447.43394,376.9991397,1.0};
    static double e[5]={ 1.0,-0.1098628627e-02,0.2734510407e-04,
  -0.2073370639e-05,0.2093887211e-06};
    static double f[5]={ -0.1562499995e-01,0.1430488765e-03,
  -0.6911147651e-05,0.7621095161e-06,-0.934935152e-07};
    static double g[5]={ 1.0,0.183105e-02,-0.3516396496e-04,
  0.2457520174e-05,-0.240337019e-06};
    static double h[5]={ 0.4687499995e-01,-0.2002690873e-03,
  0.8449199096e-05,-0.88228987e-06,0.105787412e-06};
    t=fabs(x);
    if (n<0)
{
  n=-n;
}
    if (n!=1)
{
  if (t<8.0)
  {
   y=t*t;
   p=a[5];
   q=b[5];
   for (i=4; i>=0; i--)
   {
    p=p*y+a[i];
    q=q*y+b[i];
   }
            p=p/q;
  }
        else
  {
   z=8.0/t;
   y=z*z;
            p=e[4];
   q=f[4];
            for (i=3; i>=0; i--)
   {
    p=p*y+e[i];
    q=q*y+f[i];
   }
            s=t-0.785398164;
            p=p*cos(s)-z*q*sin(s);
            p=p*sqrt(0.636619772/t);
  }
}
    if (n==0)
{
  return(p);
}
    b0=p;
    if (t<8.0)
{
  y=t*t;
  p=c[5];
  q=d[5];
        for (i=4; i>=0; i--)
  {
   p=p*y+c[i];
   q=q*y+d[i];
  }
        p=x*p/q;
}
    else
{
  z=8.0/t;
  y=z*z;
        p=g[4];
  q=h[4];
        for (i=3; i>=0; i--)
  {
   p=p*y+g[i];
   q=q*y+h[i];
  }
        s=t-2.356194491;
        p=p*cos(s)-z*q*sin(s);
        p=p*x*sqrt(0.636619772/t)/t;
}
    if (n==1)
{
  return(p);
}
    b1=p;
    if (x==0.0)
{
  return(0.0);
}
    s=2.0/t;
    if (t>1.0*n)
{
  if (x<0.0)
  {
   b1=-b1;
  }
        for (i=1; i<=n-1; i++)
  {
   p=s*i*b1-b0;
   b0=b1;
   b1=p;
  }
}
    else
{
  m=(n+(int)sqrt(40.0*n))/2;
        m=2*m;
        p=0.0;
  q=0.0;
  b0=1.0;
  b1=0.0;
        for (i=m-1; i>=0; i--)
  {
   t=s*(i+1)*b0-b1;
            b1=b0;
   b0=t;
            if (fabs(b0)>1.0e+10)
   {
    b0=b0*1.0e-10;
    b1=b1*1.0e-10;
                p=p*1.0e-10;
    q=q*1.0e-10;
   }
            if ((i+2)%2==0)
   {
    q=q+b0;
   }
            if ((i+1)==n)
   {
    p=b1;
   }
  }
        q=2.0*q-b0;
  p=p/q;
}
    if ((x<0.0)&&(n%2==1))
{
  p=-p;
}
    return(p);
  }
//功能:第二类整数阶贝塞尔函数
//参数:n-阶数
//调用:ldbesl()
double lebesl(int n,double x)
{
int i;
double y,z,p,q,s,b0,b1;
extern double ldbesl();
static double a[6]={ -2.957821389e+9,7.062834065e+9,
  -5.123598036e+8,1.087988129e+7,-8.632792757e+4,
  2.284622733e+2};
static double b[6]={ 4.0076544269e+10,7.452499648e+8,
  7.189466438e+6,4.74472647e+4,2.261030244e+2,1.0};
static double c[6]={ -4.900604943e+12,1.27527439e+12,
  -5.153438139e+10,7.349264551e+8,-4.237922726e+6,
  8.511937935e+3};
static double d[7]={ 2.49958057e+13,4.244419664e+11,
  3.733650367e+9,2.245904002e+7,1.02042605e+5,
  3.549632885e+2,1.0};
static double e[5]={ 1.0,-0.1098628627e-02,
  0.2734510407e-04,-0.2073370639e-05,
  0.2093887211e-06};
static double f[5]={ -0.1562499995e-01,
  0.1430488765e-03,-0.6911147651e-05,
  0.7621095161e-06,-0.934935152e-07};
static double g[5]={ 1.0,0.183105e-02,
  -0.3516396496e-04,0.2457520174e-05,
  -0.240337019e-06};
static double h[5]={ 0.4687499995e-01,
  -0.2002690873e-03,0.8449199096e-05,
  -0.88228987e-06,0.105787412e-06};
if (n<0)
{
  n=-n;
}
if (x<0.0)
{
  x=-x;
}
if (x==0.0)
{
  return(-1.0e+70);
}
if (n!=1)
{
  if (x<8.0)
  {
   y=x*x;
   p=a[5];
   q=b[5];
   for (i=4; i>=0; i--)
   {
    p=p*y+a[i];
    q=q*y+b[i];
   }
   p=p/q+0.636619772*ldbesl(0,x)*log(x);
  }
  else
  {
   z=8.0/x;
   y=z*z;
   p=e[4];
   q=f[4];
   for (i=3; i>=0; i--)
   {
    p=p*y+e[i];
    q=q*y+f[i];
   }
   s=x-0.785398164;
   p=p*sin(s)+z*q*cos(s);
   p=p*sqrt(0.636619772/x);
  }
}
if (n==0)
{
  return(p);
}
b0=p;
if (x<8.0)
{
  y=x*x;
  p=c[5];
  q=d[6];
  for (i=4; i>=0; i--)
  {
   p=p*y+c[i];
   q=q*y+d[i+1];
  }
  q=q*y+d[0];
  p=x*p/q+0.636619772*(ldbesl(1,x)*log(x)-1.0/x);;
}
else
{
  z=8.0/x;
  y=z*z;
  p=g[4];
  q=h[4];
  for (i=3; i>=0; i--)
  {
   p=p*y+g[i];
   q=q*y+h[i];
  }
  s=x-2.356194491;
  p=p*sin(s)+z*q*cos(s);
  p=p*sqrt(0.636619772/x);
}
if (n==1)
{
  return(p);
}
b1=p;
s=2.0/x;
for (i=1; i<=n-1; i++)
{
  p=s*i*b1-b0;
  b0=b1;
  b1=p;
}
return(p);
} //功能:变型第一类整数阶贝塞尔函数
//参数:n-阶数
//调用:
double lfbesl(int n,double x)
{
int i,m;
    double t,y,p,b0,b1,q;
    static double a[7]={ 1.0,3.5156229,3.0899424,1.2067492,
  0.2659732,0.0360768,0.0045813};
    static double b[7]={ 0.5,0.87890594,0.51498869,
  0.15084934,0.02658773,0.00301532,0.00032411};
    static double c[9]={ 0.39894228,0.01328592,0.00225319,
  -0.00157565,0.00916281,-0.02057706,
  0.02635537,-0.01647633,0.00392377};
    static double d[9]={ 0.39894228,-0.03988024,-0.00362018,
  0.00163801,-0.01031555,0.02282967,
  -0.02895312,0.01787654,-0.00420059};
    if (n<0)
{
  n=-n;
}
    t=fabs(x);
    if (n!=1)
{
  if (t<3.75)
  {
   y=(x/3.75)*(x/3.75);
   p=a[6];
            for (i=5; i>=0; i--)
   {
    p=p*y+a[i];
   }
  }
        else
  {
   y=3.75/t;
   p=c[8];
            for (i=7; i>=0; i--)
   {
    p=p*y+c[i];
   }
            p=p*exp(t)/sqrt(t);
  }
}
    if (n==0)
{
  return(p);
}
    q=p;
    if (t<3.75)
{
  y=(x/3.75)*(x/3.75);
  p=b[6];
        for (i=5; i>=0; i--)
  {
   p=p*y+b[i];
  }
        p=p*t;
}
    else
{
  y=3.75/t;
  p=d[8];
        for (i=7; i>=0; i--)
  {
   p=p*y+d[i];
  }
        p=p*exp(t)/sqrt(t);
}
    if (x<0.0)
{
  p=-p;
  ]
    if (n==1)
{
  return(p);
}
    if (x==0.0)
{
  return(0.0);
}
    y=2.0/t;
t=0.0;
b1=1.0;
b0=0.0;
    m=n+(int)sqrt(40.0*n);
    m=2*m;
    for (i=m; i>0; i--)
{
  p=b0+i*y*b1;
  b0=b1; b1=p;
        if (fabs(b1)>1.0e+10)
  {
   t=t*1.0e-10;
   b0=b0*1.0e-10;
            b1=b1*1.0e-10;
  }
        if (i==n)
  {
   t=b0;
  }
}
    p=t*q/b1;
    if ((x<0.0)&&(n%2==1))
{
  p=-p;
}
    return(p);
} //功能:变型第二类整数阶贝塞尔函数
//参数:n-阶数
//调用:lfbesl();
double lgbesl(int n,double x)
{
int i;
    double y,p,b0,b1;
    static double a[7]={ -0.57721566,0.4227842,0.23069756,
  0.0348859,0.00262698,0.0001075,0.0000074};
    static double b[7]={ 1.0,0.15443144,-0.67278579,
  -0.18156897,-0.01919402,-0.00110404,-0.00004686};
    static double c[7]={ 1.25331414,-0.07832358,0.02189568,
  -0.01062446,0.00587872,-0.0025154,0.00053208};
    static double d[7]={ 1.25331414,0.23498619,-0.0365562,
  0.01504268,-0.00780353,0.00325614,-0.00068245};
    if (n<0)
{
  n=-n;
}
    if (x<0.0)
{
  x=-x;
}
    if (x==0.0)
{
  return(1.0e+70);
}
    if (n!=1)
{
  if (x<=2.0)
  {
   y=x*x/4.0;
   p=a[6];
            for (i=5; i>=0; i--)
   {
    p=p*y+a[i];
   }
            p=p-lfbesl(0,x)*log(x/2.0);
  }
        else
  {
   y=2.0/x;
   p=c[6];
            for (i=5; i>=0; i--)
   {
    p=p*y+c[i];
   }
            p=p*exp(-x)/sqrt(x);
  }
}
    if (n==0)
{
  return(p);
}
    b0=p;
    if (x<=2.0)
{
  y=x*x/4.0;
  p=b[6];
        for (i=5; i>=0; i--)
  {
   p=p*y+b[i];
  }
        p=p/x+lfbesl(1,x)*log(x/2.0);
}
    else
{
  y=2.0/x;
  p=d[6];
        for (i=5; i>=0; i--)
  {
   p=p*y+d[i];
  }
        p=p*exp(-x)/sqrt(x);
}
    if (n==1)
{
  return(p);
}
    b1=p;
    y=2.0/x;
    for (i=1; i<n; i++)
{
  p=b0+i*y*b1;
  b0=b1;
  b1=p;
}
    return(p);
} //功能:不完全贝塔(beta)函数
//描述:Bx[a,b]=Integrate[t^(a-1)*(1-t)^(b-1),{t,0,x}]/B[a,b]
//描述:B[a,b]=gamma[a]*gamma[b]/gamma[a+b]
//参数:a-参数,b-参数
//调用:lagam();
double lhbeta(double a,double b,double x)
{
double y;
    if (a<=0.0)
{
  printf("err**a<=0!");
  return(-1.0);
}
    if (b<=0.0)
{
  printf("err**b<=0!");
  return(-1.0);
}
    if ((x<0.0)||(x>1.0))
{
  printf("err**x<0 or x>1 !");
        return(1.0e+70);
}
    if ((x==0.0)||(x==1.0))
{
  y=0.0;
}
    else
{
  y=a*log(x)+b*log(1.0-x);
        y=exp(y);
        y=y*lagam(a+b)/(lagam(a)*lagam(b));
}
    if (x<(a+1.0)/(a+b+2.0))
{
  y=y*beta(a,b,x)/a;
}
    else
{
  y=1.0-y*beta(b,a,1.0-x)/b;
}
    return(y);
}

static double beta(double a,double b,double x)
{ int k;
    double d,p0,q0,p1,q1,s0,s1;
    p0=0.0; q0=1.0; p1=1.0; q1=1.0;
    for (k=1; k<=100; k++)
{
  d=(a+k)*(a+b+k)*x;
        d=-d/((a+k+k)*(a+k+k+1.0));
        p0=p1+d*p0;
  q0=q1+d*q0;
  s0=p0/q0;
        d=k*(b-k)*x;
        d=d/((a+k+k-1.0)*(a+k+k));
        p1=p0+d*p1;
  q1=q0+d*q1;
  s1=p1/q1;
        if (fabs(s1-s0)<fabs(s1)*1.0e-07)
  {
   return(s1);
  }
}
    printf("a or b too big !");
    return(s1);
}

//功能:正态分布函数
//参数:a-均值,b-方差
//调用:lcerf(),lagam(),lbgam();
double ligas(double a,double d,double x)
{
double y;
    if (d<=0.0)
{
  d=1.0e-10;
}
    y=0.5+0.5*lcerf((x-a)/(sqrt(2.0)*d));
    return(y);
} //功能:t-分布函数
//参数:n-自由度
//调用:lhbeta(),lagam();
double ljstd(double t,int n)
{
double y;
    if (t<0.0)
{
  t=-t;
}
    y=1.0-lhbeta(n/2.0,0.5,n/(n+t*t));
    return(y);
} //功能:X^2-分布函数
//参数:n-自由度
//调用:lbgam(),lagam();
double lkchi(double x,int n)
{
double y;
    if (x<0.0)
{
  x=-x;
}
    y=lbgam(n/2.0,x/2.0);
    return(y);
} //功能:F-分布函数
//参数:n1-自由度,n2-自由度
//调用:lhbeta(),lagam();
double llf(double f,int n1,int n2)
{
double y;
    if (f<0.0)
{
  f=-f;
}
    y=lhbeta(n2/2.0,n1/2.0,n2/(n2+n1*f));
    return(y);
} //功能:正弦积分
//参数:
//调用:
double lmsi(double x)
{
int n,k,jt;
    double h,t1,t2,t,s1,s2,p;
    if (x==0.0)
{
  return(0.0);
}
    h=fabs(x);
    n=1;
t1=h*(1.0+sin(x)/x)/2.0;
    s1=t1;
jt=1;
    while (jt==1)
{
  p=0.0;
        for (k=0; k<=n-1; k++)
  {
   t=(k+0.5)*h;
   p=p+sin(t)/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07)
  {
   jt=0;
  }
        else
  {
   t1=t2;
   s1=s2;
   n=n+n;
   h=0.5*h;
  }
}
    if (x<0.0)
{
  s2=-s2;
}
    return(s2);
}

//功能:正弦积分
//参数:
//调用:
double lmsi(double x)
{
int n,k,jt;
    double h,t1,t2,t,s1,s2,p;
    if (x==0.0)
{
  return(0.0);
}
    h=fabs(x);
    n=1;
t1=h*(1.0+sin(x)/x)/2.0;
    s1=t1;
jt=1;
    while (jt==1)
{
  p=0.0;
        for (k=0; k<=n-1; k++)
  {
   t=(k+0.5)*h;
   p=p+sin(t)/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07)
  {
   jt=0;
  }
        else
  {
   t1=t2;
   s1=s2;
   n=n+n;
   h=0.5*h;
  }
}
    if (x<0.0)
{
  s2=-s2;
}
    return(s2);
} //功能:余弦积分
//参数:
//调用:
double lnci(double x)
{
int n,k,jt;
    double h,t1,t2,t,s1,s2,p,r,q;
    if (x==0.0)
{
  x=1.0e-10;
}
    if (x<0.0)
{
  x=-x;
}
    r=0.57721566490153286060651;
    q=r+log(x);
h=x;
    n=1;
t1=0.5*h*(1.0-cos(x))/x;
    s1=t1;
jt=1;
    while (jt==1)
{
  p=0.0;
        for (k=0; k<=n-1; k++)
  {
   t=(k+0.5)*h;
   p=p+(1.0-cos(t))/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07)
  {
   jt=0;
  }
        else
  {
   t1=t2;
   s1=s2;
   n=n+n;
   h=0.5*h;
  }
}
    q=q-s2;
    return(q);
} //功能:指数积分
//参数:
//调用:
double loei(double x)
{
int n,k,jt;
    double h,t1,t2,t,s1,s2,p,r,q;
    if (x==0.0)
{
  x=1.0e-10;
}
    if (x<0.0)
{
  x=-x;
}
    r=0.57721566490153286060651;
    q=r+log(x);
h=x;
    n=1;
t1=0.5*h*(-1.0+(exp(-x)-1.0)/x);
    s1=t1;
jt=1;
    while (jt==1)
{
  p=0.0;
        for (k=0; k<=n-1; k++)
  {
   t=(k+0.5)*h;
   p=p+(exp(-t)-1.0)/t;
  }
        t2=(t1+h*p)/2.0;
        s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07)
  {
   jt=0;
  }
        else
  {
   t1=t2;
   s1=s2;
   n=n+n;
   h=0.5*h;
  }
}
    q=q+s2;
    return(q);
}

//功能:第一类椭圆积分
//参数:
//调用:
static double fk(double k,double y);
double lpfk(double k,double f)
{ int n;
    double pi,y,e,ff;
    if (k<0.0)
{
  k=-k;
}
    if (k>1.0)
{
  k=1.0/k;
}
    pi=3.141592653589793;
y=fabs(f);
    n=0;
    while (y>=pi)
{
  n=n+1;
  y=y-pi;
}
    e=1.0;
    if (y>=pi/2.0)
{
  n=n+1;
  e=-e;
  y=pi-y;
}
    if (n==0)
  ff=fk(k,y);
    else
{
  ff=fk(k,pi/2.0);
        ff=2.0*n*ff+e*fk(k,y);
}
    if (f<0.0)
{
  ff=-ff;
}
    return(ff);
}

static double fk(double k,double y)
{
int n,i,jt;
    double h,t,t1,t2,s1,s2,p,q;
    if ((1.0-k<1.0e-20)&&(fabs(sin(y)-1.0)<1.0e-20))
{
  return(1.0e+10);
}
    n=1;
h=y;
    q=sqrt(1.0-k*k*sin(y)*sin(y));
    if (q<=1.0e-35)
{
  q=1.0e+35;
}
    else
{
  q=1.0/q;
}
    t1=0.5*h*(1.0+q);
    s1=t1;
jt=1;
    while (jt==1)
{
  p=0.0;
        for (i=0; i<=n-1; i++)
  {
   t=(i+0.5)*h;
            q=sqrt(1.0-k*k*sin(t)*sin(t));
            if (q<=1.0e-35)
   {
    q=1.0e+35;
   }
            else
   {
    q=1.0/q;
   }
            p=p+q;
  }
        t2=(t1+h*p)/2.0;
  s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<fabs(s2)*1.0e-07)
  {
   jt=0;
  }
        else if (h<1.0e-02)
  {
   jt=0;
  }
        else
  {
   t1=t2;
   s1=s2;
   n=n+n;
   h=0.5*h;
  }
}
    return(s2);
} //功能:第二类椭圆积分
//参数:
//调用:
static double ek(double k,double y);
double lqek(double k,double f)
{
int n;
    double pi,y,e,ff;
    if (k<0.0)
  k=-k;
    if (k>1.0)
  k=1.0/k;
    pi=3.1415926;
y=fabs(f);
    n=0;
    while (y>=pi)
{
  n=n+1;
  y=y-pi;
}
    e=1.0;
    if (y>=pi/2.0)
{
  n=n+1;
  e=-e;
  y=pi-y;
}
    if (n==0)
  ff=ek(k,y);
    else
{
  ff=ek(k,pi/2.0);
        ff=2.0*n*ff+e*ek(k,y);
}
    if (f<0.0)
  ff=-ff;
    return(ff);
}

static double ek(double k,double y)
{
int n,i,jt;
    double h,t,t1,t2,s1,s2,p,q;
    n=1;
h=y;
    q=sqrt(1.0-k*k*sin(y)*sin(y));
    t1=0.5*h*(1.0+q);
    s1=t1;
jt=1;
    while (jt==1)
{
  p=0.0;
        for (i=0; i<=n-1; i++)
  {
   t=(i+0.5)*h;
            q=sqrt(1.0-k*k*sin(t)*sin(t));
            p=p+q;
  }
        t2=(t1+h*p)/2.0;
  s2=(4.0*t2-t1)/3.0;
        if (fabs(s2-s1)<1.0e-07)
   jt=0;
        else if
   (h<1.0e-03) jt=0;
        else
  {
   t1=t2;
   s1=s2;
   n=n+n;
   h=0.5*h;
  }
}
    return(s2);
}

特殊函数计算手册 作者:张善杰,金建铭 著 出版时间:2011年版 内容简介   《特殊函数计算手册》较系统地阐述了各种特殊函数的定义、数学性质、算法、数表和程序。由特定微分方程的解定义的特殊函数有正交多项式(如Chebyshev、Laguerre和Hermite多项式),Gamma函数,Legendre函数类,Bessel函数(如球Bessel、变型Bessel、Ricatti-Bessel函数等),Kelvin函数,Airy函数,Struve函数,超几何函数,抛物柱函数,椭圆柱函数和旋转椭球函数;而由特定积分定义的特殊函数有误差函数、Fresnel积分、变型Fresnel积分、余弦和正弦积分、三类完全和不完全椭圆积分、Jacobi椭圆函数,以及指数积分等。各种特殊函数计算程序。《特殊函数计算手册》可供从事物理学、力学、应用数学、大气科学,电磁场工程、航空航天工程等学科工程技术、研究人员,以及高等院校理工科本科生、研究生和教师参考。 目录 序言 第1章 Bernoulli和Euler 数 1.1 Bernoulli数 1.2 Euler数 1.3 数表 第2章 正交多项式 2.1 引言 2.2 Chebyshev 多项式 2.3 Laguerre 多项式 2.4 Hermite 多项式 2.5 数值计算 2.6 数值积分应用 2.7 数表 第3章 Gamma,Beta和Psi函数 3.1 Gamma函数 3.2 Beta 函数 3 3 Psi 函数 3.4 不完全Gamma函数 3.5 不完全Beta函数 3.6 数表 第4章 Legendre函数 4.1 引言 4.2 第一类Legendre函数 4.3 第二类Legendre函数 4.4 第一类缔合Legendre函数 4.5 第二类缔合Legendre函数 4.6 任意次的Legendre函数 4.7 数表 第5章 Bessel函数 5.1 引言 5.2 和的计算 5.3 实宗量Bessel函数和的计算 5.4 复宗量Bessel函数和的计算 5.5 任意阶、复宗量的Bessel函数和的计算 5.6 计算的正确性和精度的评估 5.7 Bessel函数的零点 5.8 Lambda函数 5.9 数表 第6章 变型Bessel函数 6.1 引言 6.2 和的计算 6.3 实宗量变型Bessel函数和的计算 6.4 复宗量变型Bessel函数和的计算 6.5 任意阶、复宗量的变型Bessel函数和的计算 6.6 复宗量Hankell函数和的计算 6.7 数表 第7章 Bessel函数的积分 7.1 Bessel函数的简单积分 7.2 变型Bessel函数的简单积分 7.3 曲线和数表 第8章 球Bessel函数 8.1 球Bessel函数 8.2 Riccati-Bessel函数 8.3 变型球Bessel函数 8.4 数表 第9章 Kelvin函数 9.1 引言 9.2 数学性质 9.3 渐近展开式 9.4 数值计算 9.5 Kelvin函数的零点 9.6 数表 第10章 Airy函数 10.1 引言 10.2 数值计算 10.3 数表 第11章 Struve函数 11.1 Struve函数 11.2 变型Struve函数 11.3 数表 第12章 超几何函数和合流超几何函数 12.1 超几何函数的定义 12.2 超几何函数的数学性质 12.3 线性变换公式 12.4 超几何函数的递推关系式 12.5 可表为超几何函数特殊函数 12.6 超几何函数数值计算 12.7 合流超几何函数的定义 12.8 合流超几何函数的数学性质 12.9 合流超几何函数的递推关系式 12.10 可表为合流超几何函数特殊函数 12.11 Whittaker函数的定义 12.12 合流超几何函数数值计算 12.12 数表 第13章 抛物柱函数 13.1 引言 13.2 抛物柱函数的定义 13.3 主要数学性质 13.4 级数展开式和渐近展开式 13.5 数值计算 13.6 数表 第14章 Mathieu函数 14.1 Mathieu函数的定义 14.2 展开式系数和特征值的确定 14.3 特征值的近似计算 14.4 时Mathieu函数的展开式 14.5 Mathieu函数的数学性质 14.6 变型Mathieu函数的定义 14.7 变型Mathieu函数的数学性质 14.8 数值计算 14.9 数表 第15章 旋转椭球波函数 15.1 旋转椭球座标系 15.2 旋转椭球座标系中波动方程的解 15.3 长旋转椭球角向和径向波函数的定义 15.4 展开式系数和特征值的确定 15.5 第二类长旋转椭球径向波函数小时的计算 15.6 扁旋转椭球角向和径向波函数的定义 15.7 第二类扁旋转椭球径向波函数小时的计算 15.8 数值计算 15.9 数表 第16章 误差函数和Fresnel积分 16.1 误差函数定义 16.2 数值计算 16.3 Gauss概率积分 16.4 Fresnel 积分引言 16.5 Fresnel 积分的幂级数和渐近展开式 16.6 Fresnel 积分的数值计算 16.7 误差函数和Fresnel 积分的零点 16.8 数表 第17章 Cosine和Sine积分 17.1 引言 17.2 幂级数展开式和渐近展开式 17.3 数值计算 17.4 数表 第18章 椭圆积分和Jacobi椭圆函数 18.1 椭圆积分简介 18.2 椭圆积分的级数展开式 18.3 椭圆积分的数值计算 18.4 Jacobi椭圆函数引言 18.5 Jacobi椭圆函数数值计算 18.6 数表 第19章 指数积分 19.1 引言 19.2 级数展开式和连分式表示式 19.3 有理分式近似式 19.4 数值计算 19.5 数表 第20章 特殊函数计算方法综述 附录 A 几个特殊微分方程的推导 A1 Helmholtz方程和分离变量 A2 圆柱座标系 A3 椭圆柱座标系 A4 抛物柱座标系 A5 旋转椭球座标系 A6 长旋转椭球座标系 A7 扁旋转椭球座标系 A8 抛物座标系 附录 B 非线性方程的求根法 B1 Newton迭代法 B2 改进的Newton迭代法 B3 弦截法 参考文献
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值