(整理)复化梯形法复化矩形法变步长梯形变步长辛普森

陕西科技大学机械教改班用C++的积分其实积分的思想就是,微分—>求和—>取极限，如果是用纯手工法那就是先对一个函数微分，再求出它的面积，在取极限，因为我们的计算速度和计算量有限，现在有了计算机这个速度很快的机器，我们可以把微分后的每个小的面积加起来，为了满足精度，我们可以加大分区，即使实现不了微分出无限小的极限情况，我们也至少可以用有限次去接近他，下面我分析了四种不同的积分方法，和一个综合通用程序。

一．积分的基本思想1、思路：微分—>求和—>取极限。

2、Newton —Leibniz 公式 ⎰-=ba a Fb F dx x f )()()( 其中，)(x F 被积函数)(x f 的原函数。

3、用计算机积分的思路在积分区间内“微分—>求和—>控制精度”。

因为计算机求和不可以取极限，也就是不可以无限次的加下去，所以要控制精度。

二．现有的理论1、一阶求积公式---梯形公式⎰=+-=b a T b f a f a b dx x f )]()([2)( 他只能精确计算被积函数为0、1次多项式时的积分。

2、二阶求积分公式——牛顿、科特斯公式 ⎰=+++-=ba Sb f a b f a f a b dx x f )]()2(4)([6)(他只能精确计算被积函数为0、1、2、3次多项式时的积分。

三．四种实现方法1.复化矩形法将积分区间[a,b]等分成n 个子区间:],[],[],[],[],[112322110n n n n x x x x x x x x x x ---、、、 则h=(b-a)/n,区间端点值k x =a+kh)hf(x ))f(x x (x I 11121=-=)()()x (22232x hf x f x I =-=............................)()()(111n ---=-=n n n n x hf x f x x I∑==ni i x hf T 1n )(源程序：#include <iostream.h>#include<math.h>double f(double x) //计算被积函数{double y;y=log(1+x)/(1+x*x); //被积函数return y;}double Tn(double a,double b,int n) //求Tn{double t=0.0;double xk; //区间端点值double t1,t2; //用来判断精度do{double h=(b-a)/n;for(int k=1;k<=n-1;k++) //每一小段的矩形叠加{t1=t;xk=a+k*h;t+=h*f(xk);t2=t;}n++; //如果精度不够就对区间再次细分，直到达到精度要求 }while(fabs(t1-t2)<=1e-7); //判断计算精度return t;}void main(){double a=0.0; //积分下线double b=2.0; //积分上限int n=1024; //把区间分为1024段cout<<Tn(a,b,n)<<endl; //输出积分结果}执行结果：2.复化梯形法方法和复化矩形法类似，只是把原来的矩形小面积变成了梯形小面积，但是精确度明显提高了，也就是说达到同样的精度需要的时间少了。

⎥⎦⎤⎢⎣⎡++=∑∑-==111)(2)()(2n k k ni i x f b f a f h I 变形一下：∑-=++=11n )()]()([21n k k x f h b f a f h T 源程序：#include <iostream.h>#include<math.h>double f(double x) //计算被积函数{double y;y=log(1+x)/(1+x*x); //被积函数return y;}double Tn(double a,double b,int n) //求Tn{double t=0.0;double xk; //区间端点值double t1,t2,h=(b-a)/n; //用来判断精度do{h=(b-a)/n;for(int k=1;k<=n-1;k++) //余项叠加，相当于每一个小梯形相加 {t1=t;xk=a+k*h;t+=f(xk);t2=t;}n++; //如果精度不够就对区间再次细分，直到达到精度要求}while(fabs(t1-t2)<=1e-7); //判断计算精度t=h*(f(a)+f(b))/2+t*h; //加上初项就是积分结果了return t;}void main(){double a=0.0; //积分下线double b=2.0; //积分上线int n=1024; //把区间分为1024段cout<<Tn(a,b,n)<<endl; //输出积分结果}执行结果：3.变步长梯形法上面我们在应用复化积分法的时候会对区间进行分割，但是在要求精度是我们也不知道事先应该分割成多少段，分的少了精度达不到，分的多了计算量太大，为了解决这个问题我在上面加了一个循环，来不断地对区间进行再次划分，直到达到精度要求，变步长积分法与我用的方法类似，只不过变步长积分法的区间划分时成倍增加的。

实现方法；由复化梯形法知道；∑-=++=11n )()]()([21n k k x f h b f a f h T 步长h=(b-a)/n现在将步长分半，即在相邻两节点[]1,+k k x x 的中点)(21121+++=k k k x x x 处增加一个求积节点，那么 []∑∑-=+-=⋅+⋅++=1021112)(2)(24)()(n k k n n k k n n n x f h x f h b f a f h T 变形一下：n n k n n k k n n n h k h a x na b h x f h T T ⋅++=-=⋅+=+-=+∑2/)()(222110212其中：源程序：#include <iostream.h>#include <math.h>double f(double x) //计算被积函数的值{double y;y=log(1+x)/(1+x*x);return y;}double t2ntn(double a,double b){int n=1;double hn=b-a; //原步长double tn=0.0;double t2n=(f(a)+f(b))*hn/2.0;while(fabs(t2n-tn)>1e-7) //判断精度{tn=t2n;t2n=0.0;for(int k=0;k<=n-1;k++) //循环叠加t2n+=f(a+hn/2.0+k*hn);t2n=tn/2.0+t2n*hn/2.0;n=n*2;hn=hn/2.0; //步长分半}return t2n;}void main(){double a=0.0;double b=2.0;cout<<t2ntn(a,b)<<endl;}执行结果：4.变步长辛普森法之前的积分斜边都是直线，如果用抛物线接近就会更准确复化辛普森求积公式h k a x na b h x f x f b f a f h T k n k k n k k n ⋅+=-=⎥⎦⎤⎢⎣⎡+++=∑∑-=+-=/)()(4)(2)()(6102/111然后就只要每次让他的积分区间加倍就行直到达到要求的精度#include <stdio.h>#include <iostream.h>#include <math.h>double a=0.0,b=2.0,e=1e-7; //积分上下线，和精度要求int n=1024;double f(double x){Double y;y=log(1+x)/(1+x*x); //被积函数return y;}float simpson(){int i;double s,sn,s2n,p,x,h;h=(b-a)/n; //步长s=f(a)+f(b);p=0;x=a+h; //积分端点for(i=1;i<=n-1;i++){if((i%2)!=0){p=p+4*f(x); //在区间中间时乘4x=x+h;}else{s=s+2*f(x); //积分端点时乘2x=x+h;}}s=s+p; //第一次求和s2n=s*h/3; //积分值do{sn=s2n;x=a+h/2; //变步长s=s-p/2;p=0;for(i=1;i<=n;i++) //变步长只需要加就行了{p=p+4*f(x);x=x+h;}s=s+p;n=n*2;h=h/2;s2n=s*h/3;}while(fabs(s2n-sn)>e); //控制精度return s2n;}void main(){cout<<simpson()<<endl;}执行结果：四．用C++写的综合程序#include <stdio.h>#include <iostream.h>#include <math.h>class Bjhs //抽象类{public:virtual double f(double x)=0; //虚计算被积函数virtual void print()=0; //输出函数virtual double Tn()=0; //虚函数};class Fhjx:public Bjhs //复化矩形法类{public:Fhjx(){a=0;b=2;e=1e-7;n=1024;} //构造函数付初值double f(double x); //计算被积函数double Tn() //求Tn{double t=0.0;double xk; //区间端点值double t1,t2; //用来判断精度do{double h=(b-a)/n;for(int k=1;k<=n-1;k++) //每一小段的矩形叠加{t1=t;xk=a+k*h;t+=h*f(xk);t2=t;}n++; //如果精度不够就对区间再次细分，直到达到精度要求}while(fabs(t1-t2)<=e); //判断计算精度return t;}void print() //输出{cout<<"用复化矩形法计算的结果"<<Tn()<<endl;}private:double a,b,e;int n;};double Fhjx::f(double x) //外联函数{double y;y=log(1+x)/(1+x*x); //被积函数return y;}class Fhtx:public Bjhs //复化梯形法{public:Fhtx(){a=0;b=2;e=1e-7;n=1024;}double f(double x); //计算被积函数double Tn() //求Tn{double t=0.0;double xk; //区间端点值double t1,t2,h=(b-a)/n; //用来判断精度do{h=(b-a)/n;for(int k=1;k<=n-1;k++) //余项叠加，相当于每一个小梯形相加{t1=t;xk=a+k*h;t+=f(xk);t2=t;}n++; //如果精度不够就对区间再次细分，直到达到精度要求}while(fabs(t1-t2)<=1e-7); //判断计算精度t=h*(f(a)+f(b))/2+t*h; //加上初项就是积分结果了return t;}void print(){cout<<"用复化梯形法计算的结果"<<Tn()<<endl;}private:double a,b,e;int n;};double Fhtx::f(double x){double y;y=log(1+x)/(1+x*x); //被积函数return y;}class Bbctx:public Bjhs //变步长梯形法{public:Bbctx(){a=0;b=2;e=1e-7;n=1;tn=0;}double f(double x); //计算被积函数的值double Tn(){double hn=b-a; //原步长double t2n=(f(a)+f(b))*hn/2.0;while(fabs(t2n-tn)>e) //判断精度{tn=t2n;t2n=0.0;for(int k=0;k<=n-1;k++) //循环叠加t2n+=f(a+hn/2.0+k*hn);t2n=tn/2.0+t2n*hn/2.0;n=n*2;hn=hn/2.0; //步长分半}return t2n;}void print(){cout<<"用变步长梯形法计算的结果"<<Tn()<<endl;}private:double a,b,e;double tn;int n;};double Bbctx::f(double x){double y;y=log(1+x)/(1+x*x); //被积函数return y;}class Bbcxps:public Bjhs //变步长辛普森法{public:Bbcxps(){n=1024;a=0;b=2;e=1e-7;} //积分上下线，和精度要求double f(double x);double Tn(){int i;double s,sn,s2n,p,x,h;h=(b-a)/n; //步长s=f(a)+f(b);p=0;x=a+h; //积分端点for(i=1;i<=n-1;i++){if((i%2)!=0) //判奇偶，也就是看哪点乘几{p=p+4*f(x); //在区间中间时乘4x=x+h;}else{s=s+2*f(x); //积分端点时乘2x=x+h;}}s=s+p; //第一次求和s2n=s*h/3; //积分值do{sn=s2n;x=a+h/2; //变步长s=s-p/2;p=0;for(i=1;i<=n;i++) //变步长只需要加就行了{p=p+4*f(x);x=x+h;}s=s+p;n=n*2;h=h/2;s2n=s*h/3;}while(fabs(s2n-sn)<=e); //控制精度return s2n;}void print(){cout<<"用变步长辛普森法计算的结果"<<Tn()<<endl;}private:double a,b,e,p;int n;};double Bbcxps::f(double x){double y;y=log(1+x)/(1+x*x); //被积函数return y;}void display(Bjhs&q) //输出函数{q.Tn();q.print();};void main() //选择你想要的方法，用switch语句{char a;Fhjx b;Fhtx c;Bbctx d;Bbcxps e;cout<<"求函数y=log(1+x)/(1+x*x)的积分值"<<endl;cout<<"请输入你要选用的方法的编号"<<endl;cout<<"A.复化矩形法 B.复化梯形法"<<endl;cout<<"C.变步长梯形法 D.变步长辛普森法"<<endl;cin>>a;switch(a){case'A': display(b);break;case'B': display(c);break;case'C': display(d);break;case'D': display(e);break;default:cout<<"输入代码错误"<<endl;break;}}执行结果：。