数值积分计算 梯形，Simpson，Romberg

#include <vector>
#include <iostream>
#include <cmath>
using namespace std;

#define e 2.718281828459045

 double f(double x, int id){
    if(id==1) return sqrt(4.0-sin(x)*sin(x));
    if(id==2)
    {
        if(x==0) return 1;
        return sin(x)/x;
    }
    if(id==3) return pow(e,x)/(4+x*x);
    if(id==4) return log(1+x)/(1+x*x);
    cout<<"f(x)error"<<endl;
 }
double TiXing(double a,double b,int n,int id){     //分成n等分，每份用梯形面积近似
    double ans=0;
    double h=(b-a)/n;
    ans+=0.5*(f(a,id)+f(b,id));
    for(int i=1;i<=n-1;++i){
        ans+=f(a+i*h,id);
    }
    return h*ans;
}

double Simpson(double a, double b,int n,int id){   //分成2m等分，每2等分Simpson
    double ans=0,h=(b-a)/n;
    if(n%2!=0) cout<<"n必须是偶数"<<endl;
    int m=n/2;
    ans+=f(a,id)+f(b,id);
    double temp=0;
    for(int i=0;i<=m-1;++i){
        temp+=f(a+(2*i+1)*h,id);
    }
    ans+=4*temp;
    temp=0;
    for(int i=1;i<=m-1;++i){
        temp+=f(a+(2*i)*h,id);
    }
    ans+=2*temp;
    return ans*h/3;
}
int Romberg_k;
double Romberg(double a,double b,double eps,int id){
    int k;
    double h=b-a;
    vector<vector<double>>R(32,vector<double>(32,0));
    R[1][1]=(f(a,id)+f(b,id))*h/2;      //j=1时梯形公式            
    for(k=2;k<100;++k){
        R[k][1]=R[k-1][1];
        for(int i=1;i<=pow(2,k-2)+1e-5;++i)       //k=2时2等分，将[a,b] pow(2,k-2)等分
        {
            R[k][1]+=h/pow(2,k-2) * f(a+(2*i-1)*h/pow(2,k-1),id);   //hn=h/pow(2,k-2)
        }
        R[k][1]/=2;
        for(int j=2;j<=k;++j)
        {
            R[k][j]=R[k][j-1]+(R[k][j-1]-R[k-1][j-1])/(pow(4,j-1)-1);
        }
        if(fabs(R[k][k]-R[k-1][k-1])<eps) break;
    }
    Romberg_k=k;
    return R[k][k];
}

signed main()
{
    printf("%.10f\n",Romberg(0,0.25,pow(0.1,5),1));
    cout<<Romberg_k<<endl;
    system("pause");
    return 0;
}