龍貝格演算法c語言實現

Ⅰ 龍貝格求積，c語言代碼

#include<iostream.h>
#include<math.h>
# define Precision 0.000001//積分精度要求
# define e 2.71828183
#define MAXRepeat 10 //最大允許重復

double function(double x)//被積函數
{
double s;
s=pow(e,x)*cos(x);
return s;
}
double Romberg(double a,double b,double f(double x))
{
int m,n,k;
double y[MAXRepeat],h,ep,p,xk,s,q;
h=b-a;
y[0]=h*(f(a)+f(b))/2.0;//計算T`1`(h)=1/2(b-a)(f(a)+f(b));
m=1;
n=1;
ep=Precision+1;
while((ep>=Precision)&&(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=(y[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)
q=(s*p-y[k-1])/(s-1.0);//[pow(4,m)T`m`(h/2)-T`m`(h)]/[pow(4,m)-1],2m階牛頓柯斯特公式，即龍貝格公式
y[k-1]=p;
p=q;
}
ep=fabs(q-y[m-1]);//前後兩步計算結果比較求精度
m=m+1;
y[m-1]=q;
n=n+n; // 2 4 8 16
h=h/2.0;//二倍分割區間
}
return q;
}

main()
{
double a,b,Result;
cout<<"請輸入積分下限:"<<endl;
cin>>a;
cout<<"請輸入積分上限:"<<endl;
cin>>b;
Result=Romberg( a, b, function);
cout<<"龍貝格積分結果："<<Result<<endl;
return 0;
}

本文來自CSDN博客，轉載請標明出處：http://blog.csdn.net/liuhongyi666/archive/2008/12/23/3589002.aspx