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

dormouse_buaa

//Monte Carlo法计算定积分
//a-积分下限
//b-积分上限，要求b>a
//fmtclf-积分函数
//调用函数 mrnd1()
double mrnd1(double *r);
double fmtcl(double a,double b,double (*fmtclf)(double))
{
int m;
    double r,d,x,s;
    r=1.0;
s=0.0;
d=10000.0;
    for (m=0; m<=9999; m++)
{
  x=a+(b-a)*mrnd1(&r);
        s=s+fmtclf(x)/d;
}
    s=s*(b-a);
    return(s);
}
double mrnd1(double *r)
{
int m;
    double s,u,v,p;
    s=65536.0;
u=2053.0;
v=13849.0;
    m=(int)(*r/s);
*r=*r-m*s;
    *r=u*(*r)+v;
m=(int)(*r/s);
    *r=*r-m*s;
p=*r/s;
    return(p);
}

dormouse_buaa

<>//有理分式法计算二重积分
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bjpqf-积分函数
//bjpqfs-计算积分上下限y1(x),y0(x).要求y1(x)&gt;y0(x),
//其中参数y[0]返回y0(x)的值，参数y[1]返回y1(x)的值
//调用函数 ppg1()
static double pqg1(double x,double eps,double (*bjpqf)(double,double),void (*bjpqfs)(double,double*));
double bjpq(double a,double b,double eps,double (*bjpqf)(double,double),void (*bjpqfs)(double,double*))
{
int m,n,k,l,j;
double bb[10],h[10],hh,s1,s2,t1,t2,x,g,s0,ep,s;
m=1;
n=1;
hh=b-a; h[0]=hh;
s1=pqg1(a,eps,bjpqf,bjpqfs);
s2=pqg1(b,eps,bjpqf,bjpqfs);
t1=hh*(s1+s2)/2.0;
s0=t1;
bb[0]=t1;
ep=1.0+eps;
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
{
  t2=0.5*t1;
  for (k=0;k&lt;=n-1;k++)
  {
   x=a+(k+0.5)*hh;
   s1=pqg1(x,eps,bjpqf,bjpqfs);
   t2=t2+0.5*s1*hh;
  }
  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;
  }
  s=bb[m-1];
  for (j=m;j&gt;=2;j--)
  {
   s=bb[j-2]-h[j-2]/s;
  }
  ep=fabs(s-s0)/(1.0+fabs(s));
  n=n+n;
  t1=t2;
  s0=s;
  hh=hh/2.0;
}
return(s);
}</P><>static double pqg1(double x,double eps,double (*bjpqf)(double,double),void (*bjpqfs)(double,double*))
{
int m,n,k,l,j;
double b[10],h[10],y[2],hh,t1,t2,s0,yy,g,ep,s;
m=1; n=1;
bjpqfs(x,y);
hh=y[1]-y[0];
h[0]=hh;
t1=0.5*hh*(bjpqf(x,y[0])+bjpqf(x,y[1]));
s0=t1;
b[0]=t1;
ep=1.0+eps;
while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
{
  t2=0.5*t1;
  for (k=0;k&lt;=n-1;k++)
  {
   yy=y[0]+(k+0.5)*hh;
   t2=t2+0.5*hh*bjpqf(x,yy);
  }
  m=m+1;
  h[m-1]=h[m-2]/2.0;
  g=t2;
  l=0;
  for (j=2;j&lt;=m;j++)
  {
   s=g-b[j-2];
   if (l==0)
    if (fabs(s)+1.0==1.0)
    {
     l=1;
    }
    else
    {
     g=(h[m-1]-h[j-2])/s;
    }
  }
  b[m-1]=g;
  if (l!=0)
  {
   b[m-1]=1.0e+35;
  }
  s=b[m-1];
  for (j=m;j&gt;=2;j--)
  {
   s=b[j-2]-h[j-2]/s;
  }
  ep=fabs(s-s0)/(1.0+fabs(s));
  n=n+n;
  t1=t2;
  s0=s;
  hh=0.5*hh;
}
return(s);
}
</P>

dormouse_buaa

<>//变步长simpson求积法计算二重积分
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bjpqf-积分函数
//bjpqfs-计算积分上下限y1(x),y0(x).要求y1(x)&gt;y0(x),
//其中参数y[0]返回y0(x)的值，参数y[1]返回y1(x)的值
//调用函数 simp1()
static double simp1(double x,double eps,double (*bisimpf)(double,double),void (*bisimpfs)(double,double*));
double bisimp(double a,double b,double eps,double (*bisimpf)(double,double),void (*bisimpfs)(double,double*))
{
    int n,j;
    double h,d,s1,s2,t1,x,t2,g,s,s0,ep;
    n=1;
h=0.5*(b-a);
    d=fabs((b-a)*1.0e-06);
    s1=simp1(a,eps,bisimpf,bisimpfs);
s2=simp1(b,eps,bisimpf,bisimpfs);
    t1=h*(s1+s2);
    s0=1.0e+35;
ep=1.0+eps;
    while (((ep&gt;=eps)&amp;&amp;(fabs(h)&gt;d))||(n&lt;16))
{
  x=a-h;
  t2=0.5*t1;
        for (j=1;j&lt;=n;j++)
  {
   x=x+2.0*h;
            g=simp1(x,eps,bisimpf,bisimpfs);
            t2=t2+h*g;
  }
        s=(4.0*t2-t1)/3.0;
        ep=fabs(s-s0)/(1.0+fabs(s));
        n=n+n;
  s0=s;
  t1=t2;
  h=h*0.5;
}
    return(s);
}</P><>static double simp1(double x,double eps,double (*bisimpf)(double,double),void (*bisimpfs)(double,double*))
{
    int n,i;
    double y[2],h,d,t1,yy,t2,g,ep,g0;
    n=1;
    bisimpfs(x,y);
    h=0.5*(y[1]-y[0]);
    d=fabs(h*2.0e-06);
    t1=h*(bisimpf(x,y[0])+bisimpf(x,y[1]));
    ep=1.0+eps;
g0=1.0e+35;
    while (((ep&gt;=eps)&amp;&amp;(fabs(h)&gt;d))||(n&lt;16))
{ yy=y[0]-h;
t2=0.5*t1;
for (i=1;i&lt;=n;i++)
{ yy=yy+2.0*h;
t2=t2+h*bisimpf(x,yy);
}
g=(4.0*t2-t1)/3.0;
ep=fabs(g-g0)/(1.0+fabs(g));
n=n+n;
g0=g;
t1=t2;
h=0.5*h;
}
    return(g);
}</P>

dormouse_buaa

//////////////////////////////////////////////////////////////
//变步长梯形求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//func-积分函数
double bats(double a,double b,double eps,double (*func)(double));
//////////////////////////////////////////////////////////////
//变步长simplson求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//func-积分函数
double bcsimp(double a,double b,double eps,double (*func)(double));
//////////////////////////////////////////////////////////////
//自适应梯形求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//d-对积分区间进行分割的最小步长，当子区间的宽度小于d时，
//  即使没有满足精度要求，也不再往下进行分割
//bbptsf-积分函数
//调用函数 ppp()
double bbpts(double a,double b,double eps,double d,double (*bbptsf)(double);
//////////////////////////////////////////////////////////////
//龙格贝求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bdrombf-积分函数
double bdromb(double a,double b,double eps,double (*bdrombf)(double));
//////////////////////////////////////////////////////////////
//有理分式法求一维积分
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bhpqf-积分函数
double bhpq(double a,double b,double eps,double (*bhpqf)(double));
//////////////////////////////////////////////////////////////
//勒让德-高斯求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//begausf-积分函数
double begaus(double a,double b,double eps,double (*begausf)(double));
//////////////////////////////////////////////////////////////
//切比雪夫求积法
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bfcbsvf-积分函数
double bfcbsv(double a,double b,double eps,double (*bfcbsvf)(double));
//////////////////////////////////////////////////////////////
//高振荡函数求积法
//用分部积分法计算高振荡函数的积分
//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[]);
//////////////////////////////////////////////////////////////
//拉盖尔-高斯求积法计算半无限区间（0，+inf）上的积分
//flagsf-积分函数
double flags(double (*flagsf)(double));
//////////////////////////////////////////////////////////////
//Hermite-Gauss 求积公式计算无限区间（-inf,inf）上的积分
//fhmgsf-积分函数
double fhmgs(double (*fhmgsf)(double));
//////////////////////////////////////////////////////////////
//Monte Carlo法计算定积分
//a-积分下限
//b-积分上限，要求b&gt;a
//fmtclf-积分函数
//调用函数 mrnd1()
double fmtcl(double a,double b,double (*fmtclf)(double));
//////////////////////////////////////////////////////////////
//有理分式法计算二重积分
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bjpqf-积分函数
//bjpqfs-计算积分上下限y1(x),y0(x).要求y1(x)&gt;y0(x),
//其中参数y[0]返回y0(x)的值，参数y[1]返回y1(x)的值
//调用函数 ppg1()
double bjpq(double a,double b,double eps,double (*bjpqf)(double,double),void (*bjpqfs)(double,double*));
//////////////////////////////////////////////////////////////
//变步长simpson求积法计算二重积分
//a-积分下限
//b-积分上限，要求b&gt;a
//eps-积分精度要求
//bjpqf-积分函数
//bjpqfs-计算积分上下限y1(x),y0(x).要求y1(x)&gt;y0(x),
//其中参数y[0]返回y0(x)的值，参数y[1]返回y1(x)的值
//调用函数 simp1()
double bisimp(double a,double b,double eps,double (*bisimpf)(double,double),void (*bisimpfs)(double,double*));
//////////////////////////////////////////////////////////////

cherrycl

<p>我把这个积分放到C++里了，编译没有问题，但是link出错，不知是什么问题啊，小妹是菜鸟，刚刚起步，希望大侠解救啊！程序代码是</p><p>#include &lt;math.h&gt;<br/>double bats(double a,double b,double eps,double (*func)(double))<br/>&nbsp;{ <br/>&nbsp;int n,k; <br/>&nbsp;double fa,fb,h,t1,p,s,x,t; <br/>&nbsp;fa=func(a); <br/>&nbsp;fb=func(b); <br/>&nbsp;n=1; h=b-a;<br/>&nbsp;t1=h*(fa+fb)/2.0; <br/>&nbsp;p=eps+1.0; <br/>&nbsp;while (p&gt;=eps) <br/>&nbsp;{ <br/>&nbsp;&nbsp;s=0.0; <br/>&nbsp;&nbsp;for (k=0;k&lt;=n-1;k++) <br/>&nbsp;&nbsp;{ <br/>&nbsp;&nbsp;&nbsp;x=a+(k+0.5)*h; s=s+func(x); <br/>} t=(t1+h*s)/2.0; p=fabs(t1-t); <br/>&nbsp;&nbsp;t1=t; n=n+n; h=h/2.0; <br/>} <br/>&nbsp;return(t); <br/>} </p><p>连接错误显示是</p><p>LIBCD.lib(crt0.obj) : error LNK2001: unresolved external symbol _main<br/>Debug/Cpp9.exe : fatal error LNK1120: 1 unresolved externals<br/>执行 link.exe 时出错.</p><p>谢谢大家帮助啊！</p>

axiaoqiang

厉害