目录
一、 引言…………………………………… 3
二、 算法描述……………………………… 3
三、 试验及分析…………………………… 7
四、 自己的特色…………………………… 10
五、 总结…………………………………… 13
六、 程序清单……………………………… 13
一、引言
在实际应用中,常常会遇到积分制的计算,在数学分析中，我们擦汗那个擦很难过用微积分基本定理Newton-Leibniz公式来求解积分，即 ，其中F(x)是被积函数的原函数，a,b是积分区间的两个端点值，然而，随着学习的深入，发现Newton-Leibniz公式有很大的局限性:该公式要求f(x)在区间[a,b]上可积，且有原函数F(x)存在，然而对大多数f(x)而言找原函数很难，即使找到也不能用初等函数表示，如，，；再者就是某些原函数过于复杂，如，等；此外当f(x)是由测量或计算得到的数据表表示的时候，也不能使用该公式求积分。因此，为克服上述积分过程中存在的缺点，定积分计算的数值求解能弥补上述不足，并能带来满意的结果。
关键字：Newton-Leibniz 局限性 数值积分
二、算法描述
1.梯形(trapezia)公式
如果用一个梯形来近似代替每个子区间的面积，梯形的四个顶点分别位于，，与。梯形的面积为
对于整个[a,b]区间，积分值由所有窄带的和给出
2.中值发公式
如果用一个矩形来近似代替每个子区间的面积，矩形的高应该用其中点，一间距为宽可以计算出矩形面积
然后同理梯形的计算方法可以得出积分。
3. Simpson公式
当积分的上下限相对于展开的中心点对称时，积分泰勒展开式中含有f(x)的奇数阶导数的项都将等于零。利用这一性质，我们可以在相邻的两个子区间内对面积作泰勒级数展开，可得
积分的精确度可以达到。利用中心差分，可将f(x)的二阶导数近似表为
将此结果代入前述中心查分表式，则子区间[xi-1,xi+1]内的积分近似为
对于整个[a,b]区间，结果为
假设区间的数目N为偶数，则奇数格点的贡献是偶数格点贡献的2倍。这种权重上的差异来自我们为修正基本方法所得的一级结果所引入的f(x)二阶导数的贡献。而两端点和的权重则仅为偶数点的一半。Simpson法则的另一优点在于其自然引出了一种算法，即通过迭代使积分达到所需要的精确度。
误差估计---原则上，我们可以通过估计所省略的高阶项大小得出数值计算中的误差。其前提是知道高阶导数（如）的平均值），这样做并不容易。比较相继的两次迭代结果是较为简便的一种方法。例如，如果两次计算值I[a,b]之差小于预设的允许值，则我们认为计算已经达到了所需的精确度。对于Simpson法则，该算法可通过如下方法实施。假设我们在x=a和x=b之间对函数f(x)进行积分，且要求达到精确度。各个节点上的函数值对积分的贡献可分为三个部分，即端点区，奇数点与偶数点的贡献，
其中，。可用，及以如下形式写出积分
由此可见，对于一个新的N，我们可利用上一次循环中已有的，及的值。为便于进行迭代计算，应从N=开始，其中应为某个较小而合理的区间数目，例如取=6。这是迭代计算的第一步。在下一轮迭代中，将区间的数目翻倍，而令步长h减半，即令N==2。以这种方式改变子区间数目的好处在于，对于一个新的N，来自两个端点的贡献（存储为）是不变的，而偶数点与上一次迭代计算过的偶数点完全相同。也就是说，在新的迭代中的的值就是上一次循环的与之和。所需计算的仅仅是函数在（新的）奇数点上的值。从而我们得到结果
其中。利用上面等式进行变步长的循环计算即。如果计算值与前一循环所得结果之差小于ε，则可认为计算结果已收敛；否则即再次令子区间数目加倍，使N=N3=2N2。然后依此类推，直到计算结果收敛或者达到最大迭代次数
4. Gruss Qusdrature及Gruss Qusdrature-three point公式
在高斯求积公式（1-1）中，若取权函数，积分区间得
（2-1）
称之为高斯-勒让德求积公式.
对任意求积区间，通过变换
可化为区间，这时
取，，则得公式
.
公式(2-1)中求积系数
求积公式的高斯点就是勒让德多项式的零点.
定理设，求积公式（2-1）的误差为
(2-2)
利用勒让德多项式的一个性质
可得高斯-勒让德求积系数为
， （2-3）
按（2-2）式，可推得余项为
若取
的零点为节点，则
从而一点高斯-勒让德求积公式（中矩形求积公式）为
其余项为
若取
的两个零点为节点，则
从而二点高斯-勒让德求积公式为
其余项
同理，三点高斯-勒让德求积公式为
其余项
一般地，高斯-勒让德求积公式的节点可以通过勒让德多项式的零点确定，而系数通过（2-3）式确定.
三、实验及分析
1、trapezium rule
2、Mid-point rule
3、Simpson's rule
4、Gruss Qusdrature
5、Gruss Qusdrature-three point
6、下面是对五种方式在一起的一个直观的比较
四、自己的特色
特色1
代码的简洁化，代码用一个双重循环，两个主要变量实现了计算过程，在空间和时间上做到了最优化。
for (int i=0;i<20;i++)
{
area1=0;
x=PI/intervals[i];
for (int j=0;j<intervals[i];j++)
{
area1=area1+x*(f(j*x)+4*f(j*x+0.5*x)+f(j*x+x))/6;
}
double e=(area1-2)/area2;
area2=area1-2;
}
2、格式化的输出让结果更加明朗的展示在大家的面前
3、采用C++中类的思想用虚基类继承的方式实现五种不同的积分方法
虚基类：
class Integration
{
public:
virtual double integra()=0;
Integration(){x1=0;x2=3.1415926;}
Integration(double t1,double t2)
{
x1=t1;
x2=t2;
}
protected:
double x1,x2;
double f(double x);//每个类只能对一个函数积分
};
继承的子类
class Trapezium_Integration:public Integration
{
public:
Trapezium_Integration(double t1,double t2)
{
x1=t1;
x2=t2;
}
double integra();
};
double Trapezium_Integration::integra()
{
int i,j=0,intervals=1;
double xx,value_integ[20],x11;
while(intervals<=524288)
{
x11=x1,value_integ[j]=0;
xx=(x2-x1)/intervals;
value_integ[j]=value_integ[j]+f(x11);//利用公式计算
for(i=1;i<intervals;i++)
{
x11=x11+xx;
value_integ[j]=value_integ[j]+f(x11)*2;
}
//x11=x11+xx;
value_integ[j]=value_integ[j]+f(x2);
value_integ[j]=value_integ[j]*xx/2;//最终积分值
j++;
intervals=intervals*2;
}
double e[20];
for(i=0;i<20;i++)
e[i]=value_integ[i]-2;
cout.precision(15);
for(j=0;j<20;j++)
cout<<value_integ[j]<<" "<<" "<<e[j]<<" "<<e[j+1]/e[j]<<endl;
return 0;
}
五、总结
此次我们组讨论取得前所未有的成就，虽然相比其他组的还有一定的差距，但是我们组已经很圆满的达到了我们组想要实现的目的。
1、大家积极的参与
这次的实习当中，在之前大家的相互的帮助下，基础相对较差的同学询问基础较好的同学，把自己的思路和动手能想相对较弱的人说一下，在我们组经行谈论的时候，每个人都拿出了自己写的代码，虽然有的人只写了三个，相比以前，大家讨论之后会去再写，局限自己的思路的情况已经是一个很到的提高
2、百花齐放、百家争鸣
这次由于大家之前都写了代码，代码的形式也是各种各样，给各个组员除了思路上的开阔以外，在实现方式上也是一种拓宽。比如陈俊桦同学在实现的过程中让我们再去深刻的体会到了虚基类的继承原理，对C++语言也是一种开拓。韩学武同学的简洁代码告诉我们如何让程序可读性最佳，彭亚琼同学在程序运行空间和时间的把我向我们展示了程序的优化等等在这次的交流中，除了思路上大家相互了解，每个人也向大家展示了不同的代码风格。
3、大家相互的交流和学习
鉴于以上两点和大家对积分的实现思路清晰的基础上，我们用了一节多一点的时间完成了小组内的交流，剩下一部分时间我们放在了两部分上：1、帮部分出现错误的同学改错，改错过程中大家都学到一些鲜为人知有容易忽略的错误。2、大家相互交流了一下写程序的心得比如如何用面向对象的思想去写代码，如何实现代码的简洁，如何定义变量，如何实现代码空间和时间的优化等等，大家受益匪浅。
六、程序清单
#include <iostream.h>
#include<cmath>
//using namespace std;
class Integration
{
public:
virtual double integra()=0;
Integration(){x1=0;x2=3.1415926;}
Integration(double t1,double t2)
{
x1=t1;
x2=t2;
}
protected:
double x1,x2;
double f(double x);//每个类只能对一个函数积分
};
double Integration::f(double x)
{
double y;
y=sin(x);
return y;
}
////////////////////////////////////////////////////////////////////////////////
////Trapezium方法
class Trapezium_Integration:public Integration
{
public:
Trapezium_Integration(double t1,double t2)
{
x1=t1;
x2=t2;
}
double integra();
};
double Trapezium_Integration::integra()
{
int i,j=0,intervals=1;
double xx,value_integ[20],x11;
while(intervals<=524288)
{
x11=x1,value_integ[j]=0;
xx=(x2-x1)/intervals;
value_integ[j]=value_integ[j]+f(x11);//利用公式计算
for(i=1;i<intervals;i++)
{
x11=x11+xx;
value_integ[j]=value_integ[j]+f(x11)*2;
}
//x11=x11+xx;
value_integ[j]=value_integ[j]+f(x2);
value_integ[j]=value_integ[j]*xx/2;//最终积分值
j++;
intervals=intervals*2;
}
double e[20];
for(i=0;i<20;i++)
e[i]=value_integ[i]-2;
cout.precision(15);
for(j=0;j<20;j++)
cout<<value_integ[j]<<" "<<" "<<e[j]<<" "<<e[j+1]/e[j]<<endl;
return 0;
}
////////////////////////////////////////////////////////////////////////////////
//Mid_point
class Mid_point_Integration:public Integration
{
public:
Mid_point_Integration(double t1,double t2)
{x1=t1;x2=t2;}
double integra();
};
double Mid_point_Integration::integra()
{
double xx,value_integ[20];
int i,j=0,intervals=1;
while(intervals<=524288)
{
value_integ[j]=0;
double x11=x1;
xx=(x2-x1)/intervals;
for(i=0;i<intervals;i++)
value_integ[j]=value_integ[j]+f(x11+(i+0.5)*xx);
value_integ[j]=value_integ[j]*xx;
//e[j]=value_integ[j]-2;
j++;
intervals=intervals*2;
}
double e[20];
for(i=0;i<20;i++)
e[i]=value_integ[i]-2;
cout.precision(15);
for(i=0;i<20;i++)
cout<<value_integ[i]<<" "<<" "<<e[i]<<" "<<e[i+1]/e[i]<<endl;
return 0;
}
//////////////////////////////////////////////////////////////////////////////
//Simpson
class Simpson_Integration:public Integration
{
public:
double integra();
Simpson_Integration(double t1,double t2)
{
x1=t1;
x2=t2;
}
};
double Simpson_Integration::integra()
{
double xx,value_integ[20];
int i,j=0,intervals=1;
while(intervals<=2048)
{
value_integ[j]=0;
double x11=x1;
xx=(x2-x1)/intervals;
value_integ[j]=value_integ[j]+f(x1);
for(i=1;i<intervals;i=i+2)
value_integ[j]=value_integ[j]+4*f(x1+i*xx);
for(i=2;i<intervals;i=i+2)
value_integ[j]=value_integ[j]+2*f(x1+i*xx);
value_integ[j]=value_integ[j]+f(x2);
value_integ[j]=value_integ[j]*xx/3;
//e[j]=value_integ[j]-2;
j++;
intervals=intervals*2;
}
double e[20];
for(i=0;i<20;i++)
e[i]=value_integ[i]-2;
cout.precision(15);
for(i=0;i<j;i++)
cout<<value_integ[i]<<" "<<" "<<e[i]<<" "<<e[i+1]/e[i]<<endl;
return 0;
}
//////////////////////////////////////////////////////////////////////////////
//Guass
class Guass_Integration:public Integration
{
public:
double integra();
Guass_Integration(double t1,double t2)
{
x1=t1;
x2=t2;
}
};
double Guass_Integration::integra()
{
double xx,value_integ[20];
int i,j=0,intervals=1;
while(intervals<=2048)
{
value_integ[j]=0;
double x11=x1;
xx=(x2-x1)/intervals;
for(i=0;i<intervals;i++)
value_integ[j]=value_integ[j]+f(x11+i*xx+(0.5-sqrt(3)/6)*xx)+f(x11+i*xx+(0.5-sqrt(3)/6)*xx);
value_integ[j]=value_integ[j]*xx/2;
j++;
intervals=intervals*2;
}
double e[20];
for(i=0;i<20;i++)
e[i]=value_integ[i]-2;
cout.precision(15);
for(i=0;i<j;i++)
cout<<value_integ[i]<<" "<<" "<<e[i]<<" "<<e[i+1]/e[i]<<endl;
return 0;
}
//////////////////////////////////////////////////////////////////////////////
//Three_Guass
class Three_Guass_Integration:public Integration
{
public:
double integra();
Three_Guass_Integration(double t1,double t2)
{
x1=t1;
x2=t2;
}
};
double Three_Guass_Integration::integra()
{
double xx,value_integ[20];
int i,j=0,intervals=1;
while(intervals<=8192)
{
value_integ[j]=0;
double x11=x1;
xx=(x2-x1)/intervals;
for(i=0;i<intervals;i=i+1)
value_integ[j]=value_integ[j]+8*f(x11+(i+0.5)*xx)+5*f(x11+(i+0.5)*xx-sqrt(3.0/5.0)*xx/2)+5*f(x11+(i+0.5)*xx+sqrt(3.0/5.0)*xx/2);
value_integ[j]=value_integ[j]*xx/18;
j++;
intervals=intervals*2;
}
double e[20];
for(i=0;i<20;i++)
e[i]=value_integ[i]-2;
cout.precision(15);
for(i=0;i<j;i++)
cout<<value_integ[i]<<" "<<" "<<e[i]<<" "<<e[i+1]/e[i]<<endl;
return 0;
}
//////////////////////////////////////////////////////////////////////////////
int main()
{
double x1,x2;
cout<<"请输入积分区间(x1<x2)"<<endl;
cout<<"x1=";cin>>x1;
cout<<"x2=";cin>>x2;
cout<<"////////////////////////////////////////////////////////////////////////////////"<<endl;
Trapezium_Integration inte1(x1,x2);
inte1.integra();
cout<<"////////////////////////////////////////////////////////////////////////////////"<<endl;
Mid_point_Integration inte2(x1,x2);
inte2.integra();
cout<<"////////////////////////////////////////////////////////////////////////////////"<<endl;
Simpson_Integration inte3(x1,x2);
inte3.integra();
cout<<"////////////////////////////////////////////////////////////////////////////////"<<endl;
Guass_Integration inte4(x1,x2);
inte4.integra();
cout<<"////////////////////////////////////////////////////////////////////////////////"<<endl;
Three_Guass_Integration inte5(x1,x2);
inte5.integra();
// inte1.integra();
return 0;
}
#include<iostream> #include<math.h> #include <iomanip> #include<math.h> #include<stdio.h> #include"Trapezium.h" //#include"Mid_Point.h" //#include"Simpson.h" //#include"Gauss_Quadrature.h" //#include"Gauss_three.h" #define Pi 3.1415926535897931 //2.0*asin(1.0) using namespace std; int main() { double x=0.0,y=Pi; printf("区间的值为 0 ~ π\n"); //printf("%lf\n",Pi); cout<<Pi<<endl; Trapezium(x,y); cout<<endl<<endl; Mid_Point(x,y); cout<<endl<<endl; Simpson(x,y); cout<<endl<<endl; Gauss_Quadrature(x,y); cout<<endl<<endl; Gauss_three(x,y); return 0; }
#include<iostream> #include<math.h> #include <iomanip> #include<math.h> #include<stdio.h> #define Pi 3.1415926535897931 using namespace std; double function(double x) { return sin(x); }
#include<iostream> #include<math.h> #include <iomanip> #include<math.h> #include<stdio.h> #include"Gauss_three.h" #define Pi 3.1415926535897931 //2.0*asin(1.0) using namespace std; //Gauss Quadrature-two point int Gauss_Quadrature(double x1,double x2) { int i,j; double s=0,x; double error1,error2=1; double Error_Ratio; int intervals[20]= {1,2,4,8,16,32,64,128,256,512}; cout<<setiosflags(ios::left)<<setw(6)<<"No."<<" " <<setw(7)<<"No.f(x)"<<" "<<setw(15)<<"Midpoint Rule" <<" "<<setw(10)<<"Error Ratio"<<endl; for(i=0; i<10; i++) { s=0; x=(x1-x2)/intervals[i]; for(j=0; j<intervals[i]; j++) s+=function(x1+(0.5-sqrt(3.0)/6)*x+j*x)+function(x1+(0.5+sqrt(3.0)/6)*x+j*x); s=s*x/2; error1=s-2; Error_Ratio=error1/error2; error2=error1; cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]*2<<" " <<setw(15)<<error1<<" "<<setw(10)<<Error_Ratio<<endl; } return 0; }
#include<iostream> #include<math.h> #include <iomanip> #include<math.h> #include<stdio.h> #include"Function.h" #define Pi 3.1415926535897931 //2.0*asin(1.0) using namespace std; //Gauss Quadrature-three point int Gauss_three(double x1,double x2) { int i,j; double s=0,x; double error1,error2=1; double Error_Ratio; int intervals[9]= {1,2,4,8,16,32,64,128,256}; cout<<setiosflags(ios::left)<<setw(6)<<"No."<<" " <<setw(7)<<"No.f(x)"<<" "<<setw(15)<<"Midpoint Rule" <<" "<<setw(10)<<"Error Ratio"<<endl; for(i=0; i<9; i++) { s=0; x=fabs(x1-x2)/intervals[i]; x=x/2; for(j=0; j<intervals[i]; j++) s=s+x/9*(8*function(2*j*x+x)+5*function(2*j*x+x-sqrt(0.6)*x) +5*function(2*j*x+x+sqrt(0.6)*x)); error1=s-2; Error_Ratio=error1/error2; error2=error1; if(i!=0) { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]*3<<" " <<setw(15)<<error1<<" "<<setw(10)<<Error_Ratio<<endl; } else { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]*3<<" " <<setw(15)<<error1<<" "<<endl; } } return 0; }
#include<iostream> #include<math.h> #include <iomanip> #include<math.h> #include<stdio.h> #include"Simpson.h" #define Pi 3.1415926535897931 //2.0*asin(1.0) using namespace std; //中值法 int Mid_Point(double x1,double x2) { int i,j; double s=0,x; double error1,error2=1; double Error_Ratio; int intervals[20]= {1,2,4,8,16,32,64,128,256,512,1024, 2048,4096,8192,16384,32768,65536,131072,262144,524288}; cout<<setiosflags(ios::left)<<setw(6)<<"No."<<" " <<setw(7)<<"No.f(x)"<<" "<<setw(15)<<"Midpoint Rule" <<" "<<setw(10)<<"Error Ratio"<<endl; for(i=0; i<20; i++) { s=0; x=fabs(x2-x1)/intervals[i]; if(x1<x2) { for(j=0; j<intervals[i]; j++) s+=function(x1+(j+0.5)*x); } else { for(j=0; j<intervals[i]; j++) s+=function(x2+(j+0.5)*x); } s=x*s; error1=s-2; Error_Ratio=error1/error2; error2=error1; if(i!=0) { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]<<" " <<setw(15)<<error1<<" "<<setw(10)<<Error_Ratio<<endl; } else { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]<<" " <<setw(15)<<error1<<" "<<endl; } } return 0; }
#include<iostream> #include<math.h> #include <iomanip> #include<math.h> #include<stdio.h> #include"Gauss_Quadrature.h" #define Pi 3.1415926535897931 //2.0*asin(1.0) using namespace std; //Simpson's rule int Simpson(double x1,double x2) { int i,j; double s=0,x; double error1,error2=1; double Error_Ratio; int intervals[20]= {1,2,4,8,16,32,64,128,256,512}; cout<<setiosflags(ios::left)<<setw(6)<<"No."<<" " <<setw(7)<<"No.f(x)"<<" "<<setw(15)<<"Simpson Rule" <<" "<<setw(10)<<"Error Ratio"<<endl; for(i=0; i<10; i++) { s=0; //n个区间 分为2n等份 x=fabs(x2-x1)/(2*intervals[i]); s=s+function(x1)+function(x2); if(x1<x2) { for(j=1; j<2*intervals[i]; j++) { if(j%2==0) { s+=2*function(x1+j*x); } else { s+=4*function(x1+j*x); } } } else { for(j=1; j<2*intervals[i]; j++) { if(j%2==0) { s+=2*function(x2+j*x); } else { s+=4*function(x2+j*x); } } } s=s*x/3; error1=s-2; Error_Ratio=error1/error2; error2=error1; if(i!=0) { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]*2+1<<" " <<setw(15)<<error1<<" "<<setw(10)<<Error_Ratio<<endl; } else { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]*2+1<<" " <<setw(15)<<error1<<" "<<endl; } } return 0; }
#include<iostream> #include<math.h> #include <iomanip> #include<math.h> #include<stdio.h> #include"Mid_Point.h" #define Pi 3.1415926535897931 using namespace std; int Trapezium(double x1,double x2)//梯形法 { int i,j; double s=0;//记和 double x; double error1,error2=1; double Error_Ratio;//收敛性 int intervals[20]= {1,2,4,8,16,32,64,128,256,512,1024, 2048,4096,8192,16384,32768,65536,131072,262144,524288};//细分 s=s+function(x1)+function(x2); cout<<setiosflags(ios::left)<<setw(6)<<"No."<<" "<< setw(7)<<"No.f(x)"<<" "<<setw(15)<<"Trapezium Rule" <<" "<<setw(10)<<"Error Ratio"<<endl; for(i=0; i<20; i++) { s=function(x1)+function(x2); x=fabs(x2-x1)/intervals[i]; if(x1<x2) { for(j=0; j<intervals[i]-1; j++) s+=2*function(x1+(j+1)*x); } else { for(j=0; j<intervals[i]-1; j++) s+=2*function(x2+(j+1)*x); } s=x*s/2; error1=s-2; Error_Ratio=error1/error2; error2=error1; if(i!=0) { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]+1<<" " <<setw(15)<<error1<<" "<<setw(10)<<Error_Ratio<<endl; } else { cout<<setiosflags(ios::left)<<setw(6)<<intervals[i] <<" "<<setw(7)<<intervals[i]+1<<" " <<setw(15)<<error1<<" "<<endl; } } return 0; }