romberg.cpp

本人的计算方法课程作业

//romberg.cpp()--ROMBERG方法计算定积分
#include<iostream>
#include<cmath>
using namespace std;
const float ERROR=0.00001;          //设定误差
float func(float x);
float Romberg(float a,float b);
void main()
{
	float a=0.3,b=0.8,R;            //给定积分下限和积分上限
	R=Romberg(a,b);
	cout<<"利用Romberg方法计算得到的定积分值为:"<<endl;
	cout<<"R="<<R<<endl;
}
//func()定义被积函数
float func(float x)
{
	float y;
	if(x==0.0) x=0.000001;
	y=(x*x*x+sin(x))/x;
	return y;
}
//Romberg()计算给定函数func()的积分，积分区间为[a,b],误差不超过ERROR
//结果返回积分值 
float Romberg(float a,float b)
{
	float h,x,S;
	float S1,S2,T1,T2,C1,C2,R1,R2,R;
	int k=2;
	h=b-a;
	T1=(func(a)+func(b))*h/2;
A2:	S=0;
	x=a+h/2;
A3:	S=S+func(x);
	x=x+h;
	if(x<b)
	    goto A3;
	T2=(T1+h*S)/2;
	S2=T2+(T2-T1)/3;
	if(k>2)
	    goto A8;
A7:	k=k+1;
	h=h/2;
	T1=T2;
	S1=S2;
	goto A2;
A8:	C2=S2+(S2-S1)/15;
	if(k>3)
	    goto A11;
A10:	C1=C2;
	goto A7;
A11:	R2=C2+(C2-C1)/63;
	if(k>4)
	    goto A14;
A13:	R1=R2;
	goto A10;
A14:	if(fabs(R2-R1)>ERROR)   //若两次积分所得值的差大于误差精度，继续计算
	     goto A13;
	R=R2;
	return R;
}