romberg

收集的C语言算法程序

Romberg求积法2007-07-16 18:41#include<iostream>
#include<math.h>
using namespace std;
#define f(x) 1/(x)                           //注意括号。被积函数

double power(double x,int y)
{
int i;
for(i=1;i<y;i++)
x*=x;                                       //注意i<y，而不是小于等于
                                             //幂次函数是什么？
return x;
}

void main()
{
int i,j,k=1,m,n=1;
double a,b,err=0.0000001;                  
double h,T[8][8],F=0;                       //T矩阵的阶如何确定才能使程序良好运行？当T的定义阶不同时，为什么结果明显不同？    
for(i=0;i<8;i++)
for(j=0;j<8;j++)
   T[i][j]=0;                             //有必要对T初始化

cout<<"Please input a and b:\n";
cin>>a;
cin>>b;

h=(b-a)/2;
T[0][0]=h*(f(a)+f(b));

while(1)
{   
     for(i=1;i<=n;i++)        F+=f(a+(2*i-1)*h);
     
     T[k][0]=T[k-1][0]/2+h*F;

for(m=1;m<=k;m++)        T[k-m][m]=(power(4,m)*T[k-m+1][m-1]-T[k-m][m-1])/(power(4,m)-1);
  
if(fabs(T[0][m]-T[0][m-1])<err)                                     //fabs函数
{cout<<"The integration is:\n"<<T[0][m]<<"\n"; break;}              //break跳出while循环
else 
{F=0;h=h/2;n=2*n;k=k+1;}
                                                            
}
}