Uva1356/uvalive3485 Bridge桥上的绳索 (数学推导,曲线长度,定积分)

题目链接:https://www.luogu.com.cn/problem/UVA1356
在一条长度为B的线段$l$上，等距截取一些端点，两个相邻端点的距离为$d$.$(d⩽D）$。每个端点处，都作了一条长度为$H$的线段，且与线段$l$垂直。在两个相邻线段之间，都存在一条全等的抛物线，且所有抛物线的长度总和为L。

给定H,D,B,L，你需要找出一个d（$d⩽D$），使得抛物线的最小值点距线段l的距离x尽可能小（$y⩾0$），输出这个x，保留两位小数。

思路:要使得桥的数目最少,那么就让每个桥的间距恰好等于D.那么桥的数目$n=(B+D-1)/D$.
那么每一条绳索的长度是$L/n$,间距是$B/n$.令$l=L/n,b=B/n.$
现在的目标是知道$l与b要求算出x$,进一步地,我们知道了$H$,只用算出$y=H-x$就可以了,那么$y$的意义就是一个二次函数$y^{\prime }=a{x}^2$,曲线长度为$l$,在$x=b/2处对应的y^{\prime }的值$.
通过微积分的知识,我们知道一个可导函数的曲线长度计算方式是这样的:
$$s={\int }_a^b\sqrt{1+{f^{\prime }(x)}^2}dx$$.我们现在知道了曲线长度$l,想要求得他对应的在横坐标等于b/2时,对应的纵坐标值,假设这个纵坐标是$y.
那么代入坐标(b/2,y)进入$f(x)=a{x}^2$,解得$a=4\ast y/b^2$
那么$f(x)=\frac{4y}{b^2}\ast x^2$
并且对$f(x)求导,f^{\prime }(x)=\frac{8y}{b^2}\ast x$
$f^{\prime }(x{)}^2=\frac{64y^2}{b^4}\ast x^2$
代入我们已知的值进入上面的式子,由于二次函数在$[(-b/2),(b/2)]的积分区域是对称的,可以转化为两倍的[0,b/2]对应积分区域的值,也就是$
$$l=2\ast {\int }_0^{b/2}\sqrt{1+f^{\prime }(x{)}^2}dx$$
$$l=2\ast {\int }_0^{b/2}\sqrt{1+\frac{16y^2}{b^4}\ast x^2}dx$$
我们似乎陷入迷茫了,这该怎么办,我们现在知道的是$l$与求出$l$的公式,我们的目标是求得一个y值,尽量满足上述式子.显然当$y$的值越大时,弧长是越大的,我们考虑用二分$y$的值来解决这题.
显然$0<=y<=H.$
如何计算上述定积分呢?使用自适应辛普森法即可.
注意n的计算,不是向下取整的,被坑了,因为是你每个桥的间隔不得大于D,当桥的数目向下取整时,实际上每个桥的距离要大于D.所以是向上取整的.
AC代码:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
double b;double y;
double F(double x){
	return sqrt(1+ ((64*y*y)/((b*b*b*b)))*x*x );
}
double simpson(double len,double fL,double fM,double fR){
	return (len)*(fL+fR+4*fM)/6;
}
//l,r,mid都是实际的横坐标值, fL , fM ,fR 是对应的函数值 
//A是代表当前区间的辛普森值 
double asr(double l,double r,double eps,double fL,double fM,double fR,double A){
	double mid =(l+r)/2;double len = (r-l);
	double L_mid = F(l+len/4);double R_mid = F(r-len/4);
	double L=simpson(len/2,fL,L_mid,fM);double R = simpson(len/2,fM,R_mid,fR);
	//计算出左右区间的辛普森值
	if(fabs(L+R-A)<=15*eps) return (L+R)+(L+R-A)/15.0;
	else return asr(l,mid,eps/2,fL,L_mid,fM,L) + asr(mid,r,eps/2,fM,R_mid,fR,R); 
}
//使用下面的这个自适应辛普森来开始进行程序 
double asr(double l,double r,double eps){
	double L=F(l);double R = F(r);double MID = F((l+r)/2);
	double A = simpson(r-l,L,MID,R);
	return asr(l,r,eps,L,MID,R,A);
}
bool check(double x,double l){
	double eps = 1e-4;
	y = x;//y是全局变量,把y修改为x 
	double val = 2*asr(0,b/2,eps);
	return val < l;
}
int main(){
	int T;cin>>T;
	int kase=0;
	while(T--){
		++kase;
		int D,H,B,Len;
		scanf("%d %d %d %d",&D,&H,&B,&Len);
		int n = (B+D-1)/D;double l = Len*1.0/n;b = B*1.0/n;
		double L = 0, R = H;//准备二分y的值
		const double eps = 1e-4;
		while(fabs(R-L)>eps){
			double mid = (R+L)/2;
			if(check(mid,l)){
				L = mid;
			}
			else R = mid;
		}
		if(kase > 1) printf("\n"); 
		printf("Case %d:\n%.2lf\n",kase,H-L);
	}
return 0;}