romberg.cpp

常用的数值分析算法

#include<iostream.h>
#include<Math.h>
#include<iomanip.h>
#define a 0
#define b 1
#define pi 3.14159265358979323846
double c[5][5];

double T0k(double A,double B,double H,int N,double c[5][5])
{       double tres=0.0,txunhuan=0.0,tlin=0.0;      
	            if(N==0)
				{             
	            c[0][0]=(H*(cos((pi*A)/2)+cos((pi*B)/2)))/2;
				return c[0][0];
				}
             else
			 {
     for(int i=1;i<=pow(2,N-1);i++)           
	 {
		 txunhuan=txunhuan+cos((pi*(A+((2*i-1)*H)/pow(2,N)))/2);
     }   
         tlin=(c[N-1][0]+(H*txunhuan)/pow(2,N-1))/2;
		 c[N][0]=tlin;
		for(int j=1;j<=N;j++)
	{      if(j==1)
		{     
			    tres=(4*c[N][0]-c[N-1][0])/3;
                c[N][j]=tres;
		 }
		      else
			  {  
				  tres=(pow(4,j)*c[N][j-1]-c[N-1][j-1])/(pow(4,j)-1);
		          c[N][j]=tres;
			  }
		 }
             }
         return c[N][N];
			 
}



void main()
{ double h,comp=1.0;
  double result1=0.0;
  double result2=0.0;
  int k=0;
  h=b-a;
 for(k=1;k<5;k++)
 {  
     result2=T0k(a,b,h,k-1,c);	
	 result1=T0k(a,b,h,k,c);
	 comp=fabs(result1-result2);
	 if((comp<pow(10,-6)))
		 break;
	 else
	     continue;	
  }
for(int p=0;p<5;p++)
{
	for(int q=0;q<5;q++)
	{  
	   if(c[p][q]!=0)
       cout<<setiosflags(ios::fixed)<<c[p][q]<<" ";
	   else
		continue;
	}
    cout<<endl;
}
  cout<<"the result is:"<<result1<<" "<<pow(2,k)<<endl;


}