romberg.cpp

数值分析计算中的Romberg求积公式

#include<iostream.h>
#include <stdio.h>
#include <math.h>
#define E 1e-6
#define f(x) (cos(3.14159265358979323846*x/2))           //被积函数f(x)
double t[1000][1000];
int main()
{   
    int n,k,i,m,a,b;
	cout<<"请输入积分下限:"<<endl;
	cin>>a;
	cout<<"请输入积分上限:"<<endl;
	cin>>b;
    double h,g,p;
    h=(double)(b-a)/2;
    t[0][0]=h*(f(a)+f(b));
    k=1;
    n=1;
    do                                //Romberg算法
    {
        g=0;
        for (i=1;i<=n;i++)
            g+=f((a+((2*i-1)*h)));
        t[k][0]=(t[k-1][0]/2)+(h*g);
        for (m=1;m<=k;m++)
        {
            p=pow(4,(double)(m));
            t[k-m][m]=(p*t[k-m+1][m-1]-t[k-m][m-1])/(p-1);
        }
        m-=1;
        h/=2;
        n*=2;
        k+=1;
    }
    while (fabs(t[0][m]-t[0][m-1])>E); //自定义误差限E
	cout<<"被积函数f(x)在区间("<<a<<","<<b<<")上的近似值为:"<<endl;
    printf("%.9lf\n",t[0][m]);
   
    return 0;
}