数值计算程序大放送-数值积分
<P>选自<<徐世良数值计算程序集(C)>></P><P>每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。</P>
<P>今天都给贴出来了</P>
<P>#include "math.h"
//变步长梯形求积法
//a-积分下限
//b-积分上限，要求b>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>=eps)
{
s=0.0;
for (k=0;k<=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> //变步长simplson求积法
//a-积分下限
//b-积分上限，要求b>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>=eps)
{
p=0.0;
for (k=0;k<=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);
} //自适应梯形求积法
//a-积分下限
//b-积分上限，要求b>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,f0,f1,t0,z;
h=b-a;
t=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;
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<eps)||(h/2.0<d))
{
t=t+(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;
}
} //龙格贝求积法
//a-积分下限
//b-积分上限，要求b>a
//eps-积分精度要求
//func-积分函数
double bdromb(double a,double b,double eps,double (*bdrombf)(double))
{
int m,n,i,k;
double y,h,ep,p,x,s,q;
h=b-a;
y=h*(bdrombf(a)+bdrombf(b))/2.0;
m=1;
n=1;
ep=eps+1.0;
while ((ep>=eps)&&(m<=9))
{
p=0.0;
for (i=0;i<=n-1;i++)
{
x=a+(i+0.5)*h;
p=p+bdrombf(x);
}
p=(y+h*p)/2.0;
s=1.0;
for (k=1;k<=m;k++)
{
s=4.0*s;
q=(s*p-y)/(s-1.0);
y=p; p=q;
}
ep=fabs(q-y);
m=m+1;
y=q;
n=n+n;
h=h/2.0;
}
return(q);
} //有理分式法求一维积分
//a-积分下限
//b-积分上限，要求b>a
//eps-积分精度要求
//bhpqf-积分函数
double bhpq(double a,double b,double eps,double (*bhpqf)(double))
{
int m,n,k,l,j;
double h,bb,hh,t1,s1,ep,s,x,t2,g;
m=1;
n=1;
hh=b-a;
h=hh;
t1=hh*(bhpqf(a)+bhpqf(b))/2.0;
s1=t1;
bb=s1;
ep=1.0+eps;
while ((ep>=eps)&&(m<=9))
{
s=0.0;
for (k=0;k<=n-1;k++)
{
x=a+(k+0.5)*hh;
s=s+bhpqf(x);
}
t2=(t1+hh*s)/2.0;
m=m+1;
h=h/2.0;
g=t2;
l=0;
for (j=2;j<=m;j++)
{
s=g-bb;
if (l==0)
{
if (fabs(s)+1.0==1.0)
{
l=1;
}
else
{
g=(h-h)/s;
}
}
}
bb=g;
if (l!=0)
{
bb=1.0e+35;
}
g=bb;
for (j=m;j>=2;j--)
{
g=bb-h/g;
}
ep=fabs(g-s1);
s1=g;
t1=t2;
hh=hh/2.0;
n=n+n;
}
return(g);
} //勒让德-高斯求积法
//a-积分下限
//b-积分上限，要求b>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={-0.9061798459,-0.5384693101,0.0,
0.5384693101,0.9061798459};
static double c={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>=eps)&&(fabs(h)>s))
{
g=0.0;
for (i=1;i<=m;i++)
{
aa=a+(i-1.0)*h;
bb=a+i*h;
w=0.0;
for (j=0;j<=4;j++)
{
x=((bb-aa)*t+(bb+aa))/2.0;
w=w+begausf(x)*c;
}
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);
} //切比雪夫求积法
//a-积分下限
//b-积分上限，要求b>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={-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>=eps)&&(fabs(h)>d))
{
g=0.0;
for (i=1;i<=m;i++)
{
aa=a+(i-1.0)*h;
bb=a+i*h;
s=0.0;
for (j=0;j<=4;j++)
{
x=((bb-aa)*t+(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);
} //高振荡函数求积法
//用分部积分法计算高振荡函数的积分
//Integrate与Integrate
//a-积分下限
//b-积分上限，要求b>a
//m-被积函数中振荡函数的角频率
//n-给定积分区间两端点上f(x)的导数最高阶数+1
//fa-长度为n的数组，存放f(x)在积分区间端点x=a处各阶导数的数值
//fb-长度为n的数组，存放f(x)在积分区间端点x=b处各阶导数的数值
//s- 长度为2的数组，s(0)返回积分Integrate
// s(1)返回积分Integrate
//eps-积分精度要求
//begausf-积分函数
void bgpart(double a,double b,int m,int n,double fa[],double fb[],double s[])
{
int mm,k,j;
double sa,sb,ca,cb,sma,smb,cma,cmb;
sma=sin(m*a);
smb=sin(m*b);
cma=cos(m*a);
cmb=cos(m*b);
sa=sma; sa=cma; sa=-sma; sa=-cma;
sb=smb; sb=cmb; sb=-smb; sb=-cmb;
ca=cma; ca=-sma; ca=-cma; ca=sma;
cb=cmb; cb=-smb; cb=-cmb; cb=smb;
s=0.0;
s=0.0;
mm=1;
for (k=0;k<=n-1;k++)
{
j=k;
while (j>=4)
{
j=j-4;
}
mm=mm*m;
s=s+(fb*sb-fa*sa)/(1.0*mm);
s=s+(fb*cb-fa*ca)/(1.0*mm);
}
s=-s;
return;
} //拉盖尔-高斯求积法计算半无限区间（0，+inf）上的积分
//flagsf-积分函数
double flags(double (*flagsf)(double))
{
int i;
double x,g;
static double t={0.26355990,1.41340290,
3.59642600,7.08580990,12.64080000};
static double c={0.6790941054,1.638487956,
2.769426772,4.315944000,7.104896230};
g=0.0;
for (i=0; i<=4; i++)
{
x=t;
g=g+c*flagsf(x);
}
return(g);
} //Hermite-Gauss 求积公式计算无限区间（-inf,inf）上的积分
//fhmgsf-积分函数
double fhmgs(double (*fhmgsf)(double))
{
int i;
double x,g;
static double t={-2.02018200,-0.95857190,
0.0,0.95857190,2.02018200};
static double c={1.181469599,0.9865791417,
0.9453089237,0.9865791417,1.181469599};
g=0.0;
for (i=0; i<=4; i++)
{
x=t;
g=g+c*fhmgsf(x);
}
return(g);
}
页:
[1]
2