c语言编程龙贝格积分算法

c语言编程龙贝格积分算法
积分函数3x/(4+x^2),积分区间0,1

这个程序,我正好在学计算方法的时候写过,直接贴代码
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;
}
希望对您有所帮助!