龙贝格算法求定积分 C++ --数值分析

完整代码如下：

#include<iostream>
#include<math.h>
using namespace std;
double fun(double x)
{//return sin(x)+x*cos(x);//原函数return sin(x)/x;
}
double romberg(double a,double b,double(*f)(double),double eps){//上下限、被积函数、精度double T1,T2,S1,S2,C1,C2,R1,R2,S,x;double h=b-a;T1=h/2*((*f)(a)+(*f)(b));// cout<<T1;int k=0;while(1)
{   //k++;S=0;x=a+h/2;do{S=S+(*f)(x);x=x+h;}while(x<b);//if s>=bT2=T1/2+h/2*S;// cout<<T2;//if(fabs(T2-T1)<eps) return T2;S2=T2+1/3*(T2-T1);if(k==1){k++;h=h/2;T1=T2;S1=S2;//   cout<<S2;continue;}C2=S2+1/15*(S2-S1);if(k==2){k++;h=h/2;T1=T2;S1=S2;C1=C2;continue;}R2=C2+1/63*(C2-C1);if(k==3){k++;h=h/2;T1=T2;S1=S2;C1=C2;R1=R2;continue;}if(fabs(R2-R1)<eps){cout<<k<<endl;return R2;}k++;h=h/2;T1=T2;S1=S2;C1=C2;R1=R2;continue;
}}
int main()
{double eps,a,b;cout<<"请输入积分上下限: "<<endl;cin>>a>>b;cout<<"请输入精度："<<endl;cin>>eps;double ans=romberg(a,b,*fun,eps);cout<<ans;return 0;
}

龙贝格算法的程序框图：