龙贝格公式

/*
 * 龙贝格公式
 */

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

//double Partion(double (*fun)(double x), double a, double b, double &h, double &T1, double &T2);

double Romberg(double (*fun)(double x), double a, double b, double eplison)
{
	double h = b-a;
	double T1 = (h/2)*(fun(a) + fun(b));
	double T2 = 0;
	int k = 1;
	double S1=0, S2=0, C1=0, C2=0, R1=0, R2=0;

	double S, x;

	while(true)
	{
//		T2 = Partion(fun, a, b, h, T1, T2);

		cout << setprecision(7) << T1 << endl;
		cout << setprecision(7) << S1 << endl;
		cout << setprecision(7) << C1 << endl;
		cout << "--------------------------"<< endl;
		S = 0;
		x = a + h/2;
		while(x < b)
		{
			S += fun(x);
			x += h;
		}
		T2 = (1/2.0)*T1 + (h/2)*S;

		S2 = T2 + (1/3.0)*(T2-T1);
		if(k == 1)
		{
			++k;
			h /= 2;
			T1 = T2;
			S1 = S2;
			continue;
		}
		else // k!=1
		{
			C2 = S2 + (1/15.0)*(S2 - S1);
			if(k == 2)
			{
				C1 = C2;
				++k;
				h /= 2;
				T1 = T2;
				S1 = S2;
				continue;
			}
			else //k!=1 && k!=2
			{
				R2 = C2 + (1/63.0)*(C2 -C1);
				if(k == 3)
				{
					R1 = R2;
					C1 = C2;
					++k;
					h /= 2;
					T1 = T2;
					S1 = S2;
					continue;
				}
				else //k!=1 && k!=2 && k!=3
				{
					if(abs(R2-R1)>=eplison)
					{
						R1 = R2;
						C1 = C2;
						++k;
						h /= 2;
						T1 = T2;
						S1 = S2;
						continue;
					}
					else
					{
						break;
					}
				}
			}
		}
	}

	return R2;

}

//double Partion(double (*fun)(double x), double a, double b, double &h, double &T1, double &T2)
//{
	double h = b-a;
	double T1 = (h/2)*(fun(a) + fun(b));
//	double S, x;
	double T2 = 0;
//
//	cout << setprecision(8) << T1 << endl;
//	S = 0;
//	x = a + h/2;
//	while(x < b)
//	{
//		S += fun(x);
//		x += h;
//	}
//	T2 = (1/2.0)*T1 + (h/2)*S;
//
//	return T2;
//}

//待积分函数
double function(double x)
{
	if(x == 0)
		return 1;
	else
		return sin(x)/x;
}

int main()
{
	cout << "-----Romberg-----" << endl;
	cout << Romberg(function, 0, 1, 0.0000001) << endl;
	return 0;
}