自适应simpson积分算法.txt

这是一个很实用的求积分算法

#include <math.h>
#include <stdlib.h>
#include <fstream.h>
#include <string.h>
#include <iostream.h>
#define n 100
float f(float x)
{
	float y;
	y=x*((float)sqrt(1+x*x));
	return(y);
}
void main()
{
        int i,m;
        float APP,e[n],a[n],h[n],ga[n],gc[n],gb[n],s[n];
        float b,c,v[9],l[n],p1,p2,q,gd,ge;
 	cout<<"请输入下界A： ";
 	cin>>b;
 	cout<<"请输入上界B：  ";
 	cin>>c;
 	cout<<"请输入误差容限TOL：  ";
 	cin>>q;
 	cout<<"请输入最大基数：  ";
 	cin>>m;
 
 
 	APP=0;
 	i=1;
 	e[i]=10*q;
 	a[i]=b;
 	h[i]=(c-b)/2;
 	ga[i]=f(b);
 	gc[i]=f(a[i]+h[i]);
 	gb[i]=f(c);
 	s[i]=h[i]*(ga[i]+4*gc[i]+gb[i])/3;
 	l[i]=1;
 
 
 loop:	gd=f(a[i]+h[i]/2);
 		ge=f(a[i]+3*h[i]/2);
 		p1=h[i]*(ga[i]+4*gd+gc[i])/6;
 		p2=h[i]*(gc[i]+4*ge+gb[i])/6;
 		v[1]=a[i];
 		v[2]=ga[i];
 		v[3]=gc[i];
 		v[4]=gb[i];
 		v[5]=h[i];
 		v[6]=e[i];
 		v[7]=s[i];
 		v[8]=l[i];
 
 
 		i=i-1;//删除此级
  
 		if(fabs(p1+p2-v[7])<v[6])
 		{
 			APP=APP+p1+p2;
 		}
 		else
 		{
 			if(v[8]>m)
 			{
 				cout<<"LEVEL EXCEEDED"<<endl;
 				exit(0);
 			}
 
 			i=i+1;//右半子区间数据
 		        a[i]=v[1]+v[5];
 			ga[i]=v[3];
 			gc[i]=ge;
 			gb[i]=v[4];
 			h[i]=v[5]/2;
 			e[i]=v[6]/2;
 			s[i]=p2;
 			l[i]=v[8]+1;
 
 
 			i=i+1;//左半子区间数据
                        a[i]=v[1];
 		        ga[i]=v[2];
 			gc[i]=gd;
 			gb[i]=v[3];
 			h[i]=h[i-1];
 			e[i]=e[i-1];
 			s[i]=p1;
 			l[i]=l[i-1];
 		}
 		if(i>0) goto loop;
 		ofstream ofs("result.txt");
 
 
                ofs<<"下界A： ";
                ofs<<b<<endl;
 	        ofs<<"上界B：  ";
 	        ofs<<c<<endl;
 	        ofs<<"误差容限TOL：  ";
 	        ofs<<q<<endl;
                ofs<<"最大基数：  ";
 	        ofs<<m<<endl;
 
 
 		ofs<<"计算结果为：   ";
 		ofs<<APP;
         
 
 		cout<<"计算结果为： "<<APP<<endl;
}