romberg.cpp

龙贝格方法。又一个数值算法。

#include<stdio.h>
#include<math.h>
#include<iomanip.h>
#include<iostream.h>
double f(double x,int flag)
{
	double y;
    if (flag==1)	y=1/(1+x*x);
	if (flag==2)	y=exp(-x*x);
 	return (y);
}

double Romberg(double a,double b,double eps,int flag)
{
	int k=0;
	double h1,h2,t1,t2,s1=0,s2=0,c1=0,c2=0,r1=0,r2=0,max=10;
	double s,x;
	h1=b-a;
	t1=h1*(f(a,flag)+f(b,flag))/2;
	cout<<"k"<<setw(15)<<"T"<<setw(15)<<"S"<<setw(15)<<"C"<<setw(15)<<"R"<<endl;
	cout<<k<<setw(15)<<t1<<setw(15)<<s2<<setw(15)<<c2<<setw(15)<<r2<<endl;
	while(1)
	{
		h2=h1/2;
		x=a+h2;
		s=0;
		k++;
		do
		{
			s=s+f(x,flag);
			x=x+h1;
        }while(x<b);
		t2=(t1+h1*s)/2;
		if (k>=1)
			s2=(4*t2-t1)/3; 
		if(k>=2)
		    c2=(16*s2-s1)/15;  
		if(k>=3)
			r2=(64*c2-c1)/63; 
		if(k>=4)
		{
			max=fabs(r2-r1);
			if(max<eps)
				break;
		}
	    cout<<k<<setw(15)<<t2<<setw(15)<<s2<<setw(15)<<c2<<setw(15)<<r2<<endl;
	 	h1=h2;
		t1=t2;
		s1=s2;
		c1=c2;
		r1=r2;

	}
	return r2;
}

void main()
{
	int k=0,flag=1;
	double a,b;  
	double eps,r2;
	cout<<"第五单元 数值积分"<<endl;
	cout<<"1、用龙贝格计算积分f(x)=1/(1+x*x),x属于(0,1) "<<endl;
	cout<<"2、用龙贝格计算积分f(x)=exp(-x*x),x属于(0,0.8) "<<endl;
	cout<<"解第一题："<<endl;
	a=0;
	b=1;
    eps=0.000001;   
	r2=Romberg(a,b,eps,1);
	cout<<"利用龙贝格方法计算，其积分值为: "<<r2<<endl;
	cout<<endl<<"解第二题："<<endl;
	a=0;
	b=0.8;
	r2=Romberg(a,b,eps,2);
	cout<<"利用龙贝格方法计算，其积分值为: "<<r2<<endl;
}