龙贝格积分.cpp

数值计算的例子和参考书籍

//龙贝格积分,数值方法计算定积分
#include <cmath>
#include <vector>
/* Numerical Integration (Romberg) */
template <typename T>
T romberg(T (*f)(T), T a, T b, T eps = 1e-15)
{
	std::vector<T> R;
	int k = -1;
	T r = 0.5 * (b - a) * (f(a) + f(b));
	R.push_back(r);
	do
	{
		k += 1;
		r = 0.0;
		for(int i = 0; i < pow(2, k); i++)
		{
			r += f(a + (b - a) * (i + 0.5) / pow(2, k));
		}
		r *= (b - a) / pow(2, k + 1); 
		r += 0.5 * R[k];
		R.push_back(r);
		for(int m = 0; m <= k; m++)
		{
			R[k - m] = (pow(4, m + 1) * R[k + 1 - m] - R[k - m]) / (pow(4, m + 1) - 1);
		}
	}
	while(fabs(R[0] - R[1]) > eps);
	return R[0];
}

//例题 zju2614 (Hangzhou 2005 B)
/**
 * Problem: zju2614
 * Source : ICPC Asia Regional (Hangzhou), 2005, B Bridge
 * Author : LIU Yu, <[e-mail]>
 * Status : AC
**/

#include <cstdio>
#include <cmath>
#include <vector>

/* Numerical Integration (Romberg) */
template <typename T>
T romberg(T (*f)(T), T a, T b, T eps = 1e-15)
{
	std::vector<T> R;
	int k = -1;
	T r = 0.5 * (b - a) * (f(a) + f(b));
	R.push_back(r);
	do
	{
		k += 1;
		r = 0.0;
		for(int i = 0; i < pow(2, k); i++)
		{
			r += f(a + (b - a) * (i + 0.5) / pow(2, k));
		}
		r *= (b - a) / pow(2, k + 1);
		r += 0.5 * R[k];
		R.push_back(r);
		for(int m = 0; m <= k; m++)
		{
			R[k - m] = (pow(4, m + 1) * R[k + 1 - m] - R[k - m]) / (pow(4, m + 1) - 1);
		}
	}
	while(fabs(R[0] - R[1]) > eps);
	return R[0];
}

/* Dichotomy */
template <typename T> 
T dicho(T (*f)(T), T a, T b, T eps = 1e-15)
{
	T c = 0.5 * (a + b);
	do
	{
		if (f(c) * f(a) > 0)
		{
			a = c;
		}
		else
		{
			b = c;
		}
		c = 0.5 * (a + b); 
	}
	while(b - a > eps);
	return c;
}

template <typename T> 
T inline sqr(T a)
{
	return a * a;
}

double epsilon = 1e-6;

namespace
{
	template <double &A>
		double f(double t)
		{
			return sqrt(1.0 + 4.0 * sqr(A) * sqr(t));
		}
	template <double &A>
		double g(double width)
		{
			return 2.0 * romberg(f<A>, 0.0, 0.5 * width, epsilon);
		}
	template <int &N, double &A, double &B, double &L>
		double p(double h)
		{
			A = 4.0 * sqr(N) * h / sqr(B);
			return g<A>(B / N) - L / N;
		}
	int N;
	double A;
	double D, B, H, L;
}

void inline Bridge(void)
{
	using std::min;
	scanf("%lf %lf %lf %lf\n", &D, &H, &B, &L);
	N = (int)ceil(B / min(B, D));
	double h = dicho(p<N,A,B,L>, 0.0, H, epsilon);
	printf("%.2lf\n", H - h);
}

int main(int argc, char * argv[])
{
	int T = 0;
	scanf("%d\n", &T);
	for(int t = 0; t < T; t++)
	{
		printf("Case %d:\n", t + 1);
		Bridge();
		if(t < T - 1)
		{
			printf("\n");
		}
		fflush(stdout);
	}
	return 0;
}