具体算法详见参考文献,本文仅仅给出C++利用Romberg算法计算积分的具体程序。
romberg.h:
#ifndef ROMBERG_H_INCLUDED
#define ROMBERG_H_INCLUDED
#include <iostream>
#include <math.h>
#include"fuction.h"
using namespace std;
double Romberg(double eps,double lower_bd,double upper_bd)
{
double h=upper_bd-lower_bd,S,t1,t2,s1,s2,c1,c2;
int i,n=1,k=1;
S=fuction(lower_bd+0.5*h);
t1=0.5*h*(fuction(lower_bd)+fuction(upper_bd));
t2=t2=0.5*t1+0.5*S*h;
s2=t2+(t2-t1)/3;
c2=s2+(s2-s1)/15;
do
{
k=k+1;
h=0.5*h;
t1=t2;
s1=s2;
n=2*n;
S=0;
for(i=0; i<n; i++)
{
S=S+fuction(lower_bd+(i+0.5)*h);
}
t2=0.5*t1+0.5*S*h;
s2=t2+(t2-t1)/3;
c1=c2;
c2=s2+(s2-s1)/15;
}while (fabs(c2-c1)>eps);
return c2;
}
#endif
在fuction.h中存放具体的被积函数,下面举一例
#ifndef FUCTION_H_INCLUDED
#define FUCTION_H_INCLUDED
#include <iostream>
#include <math.h>
using namespace std;
double fuction(double x)
{
return pow(4-(pow(sin(x),2)),0.5);
}
#endif // FUCTION_H_INCLUDED
主函数中调用函数Romberg进行计算:
#include <iostream>
#include <iomanip>
#include"romberg.h"
#include"fuction.h"
using namespace std;
int main()
{
cout<<fixed<<setprecision(9)<<Romberg(1e-10,0,0.25)<<endl;
return 0;
}
经过计算,积分值为0.4987111176。
参考文献:黄云清,舒适. 2009.数值计算方法.北京:科学出版社