這個程式，我正好在學計算方法的時候寫過，直接貼程式碼

C++實現如下：

#include<iostream>

#include<cmath>

using namespace std;

const int MAXRepeat = 100; //最大允許重複

double function(double x)//被積函式,根據自己的需要手工輸入

{

double s;

s = 1.0 / (1 + x);

return s;

}

void Romberg(double a, double b, double epsion, double f(double x))

{

int m = 1;

int n = 1;

int k;

double h;

double ep;

double p;

double xk;

double s;

double q;

double T[MAXRepeat];

h = b - a;

T[0] = 0.5 * h * (f(a) + f(b));

ep = epsion + 1.0;

while ((ep >= epsion) && (m < MAXRepeat))

{

p = 0.0;

for (k = 0; k < n; k++)

{

xk = a + (k + 0.5) * h; // n-1

p = p + f(xk); //計算∑f(xk+h/2),T

} // k=0

p = (T[0] + h * p) / 2.0; //T`m`(h/2),變步長梯形求積公式

s = 1.0;

for (k = 1; k <= m; k++)

{

s = 4.0 * s; //[pow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1],2m階牛頓柯斯特公式，即龍貝格公式

q = (s * p - T[k - 1]) / (s - 1.0);

T[k-1] = p;

p = q;

}

ep = fabs(q - T[m - 1]);

m++;

T[m - 1] = q;

n++; // 2 4 8 16

h /= 2.0;

}

for (int i = 0; i < m; i++)

{

int j;

if (!(i % j))

{

cout<<T[i]<<endl;

}

else

{

cout<<T[i]<<" ";

}

j++;

}

}

int main()

{

double a;

double b;

double epsion;

cout<<"Please input the lower limit: ";

cin>>a;

cout<<"Please input the upper limit: ";

cin>>b;

cout<<"Please input the precision : ";

cin>>epsion;

Romberg( a, b, epsion, function);

return 0;

}