数值计算程序大放送-数值积分

dormouse_buaa

<>选自&lt;&lt;徐世良数值计算程序集(C)&gt;&gt;</P>
<>每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。</P>
<>今天都给贴出来了</P>
<P>#include "math.h"
//变步长梯形求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//func-积分函数
double bats(double a,double b,double eps,double (*func)(double))
{
    int n,k;
    double fa,fb,h,t1,p,s,x,t;
    fa=func(a);
    fb=func(b);
    n=1;
    h=b-a;
    t1=h*(fa+fb)/2.0;
    p=eps+1.0;
    while (p&gt;=eps)
{
  s=0.0;
        for (k=0;k&lt;=n-1;k++)
  {
   x=a+(k+0.5)*h;
            s=s+func(x);
  }
        t=(t1+h*s)/2.0;
        p=fabs(t1-t);
        t1=t;
        n=n+n;
        h=h/2.0;
}
    return(t);
}
</P>

dormouse_buaa

//变步长simplson求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//func-积分函数
double bcsimp(double a,double b,double eps,double (*func)(double))
{
int n,k;
double h,t1,t2,s1,s2,ep,p,x;
n=1;
h=b-a;
t1=h*(func(a)+func(b))/2.0;
s1=t1;
ep=eps+1.0;
while (ep&gt;=eps)
{
  p=0.0;
  for (k=0;k&lt;=n-1;k++)
  {
   x=a+(k+0.5)*h;
   p=p+func(x);
  }
  t2=(t1+h*p)/2.0;
  s2=(4.0*t2-t1)/3.0;
  ep=fabs(s2-s1);
  t1=t2;
  s1=s2;
  n=n+n;
  h=h/2.0;
}
return(s2);
}

dormouse_buaa

//自适应梯形求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//d-对积分区间进行分割的最小步长，当子区间的宽度小于d时，
//  即使没有满足精度要求，也不再往下进行分割
//bbptsf-积分函数
//调用函数 ppp()
static void ppp(double x0,double x1,double h,double f0,double f1,double t0,double eps,double d,double t[],double (*bbptsf)(double));
double bbpts(double a,double b,double eps,double d,double (*bbptsf)(double))
{
    double h,t[2],f0,f1,t0,z;
    h=b-a;
t[0]=0.0;
    f0=bbptsf(a);
f1=bbptsf(b);
    t0=h*(f0+f1)/2.0;
    ppp(a,b,h,f0,f1,t0,eps,d,t,bbptsf);
    z=t[0];
return(z);
}
static void ppp(double x0,double x1,double h,double f0,double f1,double t0,double eps,double d,double t[],double(*bbptsf)(double))
{
    double x,f,t1,t2,p,g,eps1;
    x=x0+h/2.0;
f=bbptsf(x);
    t1=h*(f0+f)/4.0;
t2=h*(f+f1)/4.0;
    p=fabs(t0-(t1+t2));
    if ((p&lt;eps)||(h/2.0&lt;d))
{
  t[0]=t[0]+(t1+t2);
  return;
}
    else
{
  g=h/2.0;
  eps1=eps/1.4;
        ppp(x0,x,g,f0,f,t1,eps1,d,t,bbptsf);
        ppp(x,x1,g,f,f1,t2,eps1,d,t,bbptsf);
        return;
}
}

dormouse_buaa

//龙格贝求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//func-积分函数
double bdromb(double a,double b,double eps,double (*bdrombf)(double))
{
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
    h=b-a;
    y[0]=h*(bdrombf(a)+bdrombf(b))/2.0;
    m=1;
n=1;
ep=eps+1.0;
    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
{
  p=0.0;
        for (i=0;i&lt;=n-1;i++)
  {
   x=a+(i+0.5)*h;
            p=p+bdrombf(x);
  }
        p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k&lt;=m;k++)
  {
   s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
  }
        ep=fabs(q-y[m-1]);
        m=m+1;
  y[m-1]=q;
  n=n+n;
  h=h/2.0;
}
    return(q);
}

dormouse_buaa

//有理分式法求一维积分
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bhpqf-积分函数
double bhpq(double a,double b,double eps,double (*bhpqf)(double))
{
    int m,n,k,l,j;
    double h[10],bb[10],hh,t1,s1,ep,s,x,t2,g;
    m=1;
n=1;
    hh=b-a;
h[0]=hh;
    t1=hh*(bhpqf(a)+bhpqf(b))/2.0;
    s1=t1;
bb[0]=s1;
ep=1.0+eps;
    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
{
  s=0.0;
        for (k=0;k&lt;=n-1;k++)
  {
   x=a+(k+0.5)*hh;
            s=s+bhpqf(x);
  }
        t2=(t1+hh*s)/2.0;
        m=m+1;
        h[m-1]=h[m-2]/2.0;
        g=t2;
        l=0;
        for (j=2;j&lt;=m;j++)
  {
   s=g-bb[j-2];
            if (l==0)
   {
    if (fabs(s)+1.0==1.0)
    {
     l=1;
    }
    else
    {
     g=(h[m-1]-h[j-2])/s;
    }
   }
  }
        bb[m-1]=g;
        if (l!=0)
  {
   bb[m-1]=1.0e+35;
  }
        g=bb[m-1];
        for (j=m;j&gt;=2;j--)
  {
   g=bb[j-2]-h[j-2]/g;
  }
        ep=fabs(g-s1);
        s1=g;
  t1=t2;
  hh=hh/2.0;
  n=n+n;
}
    return(g);
}

dormouse_buaa

//勒让德-高斯求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//begausf-积分函数
double begaus(double a,double b,double eps,double (*begausf)(double))
{
    int m,i,j;
    double s,p,ep,h,aa,bb,w,x,g;
    static double t[5]={-0.9061798459,-0.5384693101,0.0,
  0.5384693101,0.9061798459};
    static double c[5]={0.2369268851,0.4786286705,0.5688888889,
  0.4786286705,0.2369268851};
    m=1;
    h=b-a;
s=fabs(0.001*h);
    p=1.0e+35;
ep=eps+1.0;
    while ((ep&gt;=eps)&amp;&amp;(fabs(h)&gt;s))
    {
  g=0.0;
        for (i=1;i&lt;=m;i++)
        {
   aa=a+(i-1.0)*h;
   bb=a+i*h;
            w=0.0;
            for (j=0;j&lt;=4;j++)
            {
    x=((bb-aa)*t[j]+(bb+aa))/2.0;
                w=w+begausf(x)*c[j];
            }
            g=g+w;
        }
        g=g*h/2.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g;
  m=m+1;
  h=(b-a)/m;
    }
    return(g);
}

dormouse_buaa

//切比雪夫求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bfcbsvf-积分函数
double bfcbsv(double a,double b,double eps,double (*bfcbsvf)(double))
{
    int m,i,j;
    double h,d,p,ep,g,aa,bb,s,x;
    static double t[5]={-0.8324975,-0.3745414,0.0,
  0.3745414,0.8324975};
    m=1;
    h=b-a;
d=fabs(0.001*h);
    p=1.0e+35;
ep=1.0+eps;
    while ((ep&gt;=eps)&amp;&amp;(fabs(h)&gt;d))
{
  g=0.0;
        for (i=1;i&lt;=m;i++)
  {
   aa=a+(i-1.0)*h;
   bb=a+i*h;
            s=0.0;
            for (j=0;j&lt;=4;j++)
   {
    x=((bb-aa)*t[j]+(bb+aa))/2.0;
                s=s+bfcbsvf(x);
   }
            g=g+s;
  }
        g=g*h/5.0;
        ep=fabs(g-p)/(1.0+fabs(g));
        p=g;
  m=m+1;
  h=(b-a)/m;
}
    return(g);
}

dormouse_buaa

//高振荡函数求积法
//用分部积分法计算高振荡函数的积分
//Integrate[f(x)*sin(m*x)dx,{a,b}]与Integrate[f(x)*cos(m*x)dx,{a,b}]
//a-积分下限
//b-积分上限，要求b&gt;a
//m-被积函数中振荡函数的角频率
//n-给定积分区间两端点上f(x)的导数最高阶数+1
//fa-长度为n的数组，存放f(x)在积分区间端点x=a处各阶导数的数值
//fb-长度为n的数组，存放f(x)在积分区间端点x=b处各阶导数的数值
//s- 长度为2的数组，s(0)返回积分Integrate[f(x)*cos(m*x)dx,{a,b}]
//   s(1)返回积分Integrate[f(x)*sin(m*x)dx,{a,b}]
//eps-积分精度要求
//begausf-积分函数
void bgpart(double a,double b,int m,int n,double fa[],double fb[],double s[])
{
int mm,k,j;
    double sa[4],sb[4],ca[4],cb[4],sma,smb,cma,cmb;
    sma=sin(m*a);
smb=sin(m*b);
    cma=cos(m*a);
cmb=cos(m*b);
    sa[0]=sma; sa[1]=cma; sa[2]=-sma; sa[3]=-cma;
    sb[0]=smb; sb[1]=cmb; sb[2]=-smb; sb[3]=-cmb;
    ca[0]=cma; ca[1]=-sma; ca[2]=-cma; ca[3]=sma;
    cb[0]=cmb; cb[1]=-smb; cb[2]=-cmb; cb[3]=smb;
    s[0]=0.0;
s[1]=0.0;
    mm=1;
    for (k=0;k&lt;=n-1;k++)
{
  j=k;
        while (j&gt;=4)
  {
   j=j-4;
  }
        mm=mm*m;
        s[0]=s[0]+(fb[k]*sb[j]-fa[k]*sa[j])/(1.0*mm);
        s[1]=s[1]+(fb[k]*cb[j]-fa[k]*ca[j])/(1.0*mm);
}
    s[1]=-s[1];
    return;
}

dormouse_buaa

//拉盖尔-高斯求积法计算半无限区间（0，+inf）上的积分
//flagsf-积分函数
double flags(double (*flagsf)(double))
{
int i;
double x,g;
static double t[5]={0.26355990,1.41340290,
  3.59642600,7.08580990,12.64080000};
static double c[5]={0.6790941054,1.638487956,
  2.769426772,4.315944000,7.104896230};
g=0.0;
for (i=0; i&lt;=4; i++)
{
  x=t;
  g=g+c*flagsf(x);
}
return(g);
}

dormouse_buaa

//Hermite-Gauss 求积公式计算无限区间（-inf,inf）上的积分
//fhmgsf-积分函数
double fhmgs(double (*fhmgsf)(double))
{
int i;
    double x,g;
    static double t[5]={-2.02018200,-0.95857190,
  0.0,0.95857190,2.02018200};
    static double c[5]={1.181469599,0.9865791417,
  0.9453089237,0.9865791417,1.181469599};
    g=0.0;
    for (i=0; i&lt;=4; i++)
{
  x=t;
  g=g+c*fhmgsf(x);
}
    return(g);
}