LA3485二分+求解积分方程+辛普森算法计算积分

/*LA3485:
求解积分方程
关于这道题的数学模型：
给定抛物线长度L，抛物线函数f(x)=a(x-d)(x+d),求解 |a*d*d|的值，a>0,曲线积分函数lf(x)=sqrt(1+f'(x)^2)
a越大，|a*d*d|越大，L越长，所以可以二分求解

辛普森算法解函数f(x)在区间(a,b)上的积分：
模板如下：
double simpsonF(double a,double b);//返回测试值
double simpsonM(double a,double b,eps,double A);//自适应simpson递归
//答案是：simpsonM(a,b,eps,simpsonF(a,b));
double simpsonF(double a,double b)
{
double c=a+(b-a)/2;
return (F(a)+4*F(c)+F(b))*(b-a)/6;
}
double simpsonM(double a,double b,double eps,double A)
{
double c=a+(b-a)/2;
double L=simpsonF(a,c),R=simpsonF(c,b);
if (fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0;
return simpsonM(a,c,eps/2,L)+simpsonM(c,b,eps/2,R);
}
简要原理分析：细节部分自己以后需要学习！
传统求微积分的方法，是手动计算，但是很多难以写出被积函数，这里用计算机帮助计算。
任然使用先微分，再求和的方法。微分就是等分区间，去重点函数值，近似成矩形，最终将面积求和。
如果直接使用这个方法，越细分，近似程度越高，计算量越大，不细分，精度不够。以前遇到题目，用朴素算法写错了。
我们在函数曲线平缓的地方，少细分一些，在陡峭的地方（容易损失精度）多细分。折中的办法满足了时间和精度的要求。

*/


#include<iostream>
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<stdlib.h>
#include<math.h>
#include<queue>
#include<vector>
#include<map>
#define LL long long
using namespace std;
double d,l;
double F(double k,double x)//k是枚举的系数
{
    k=k/d/d;
    return sqrt(1+(2*k*(x))*(2*k*(x)));
}
double simpsonF(double a,double b,double k)
{
double c=a+(b-a)/2;
return (F(a,k)+4*F(c,k)+F(b,k))*(b-a)/6;
}
double simpsonM(double a,double b,double eps,double A,double k)
{
double c=a+(b-a)/2;
double L=simpsonF(a,c,k),R=simpsonF(c,b,k);
if (fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0;
return simpsonM(a,c,eps/2,L,k)+simpsonM(c,b,eps/2,R,k);
}
double Q(double k)
{
    double ans=simpsonM(-d,d,1e-7,simpsonF(-d,d,k),k);
    return ans-l;
}
int T,D,H,B,L;
int main()
{
//    freopen("out.txt","w",stdout);
    cin>>T;
    for(int cas=1;cas<=T;cas++)
    {
        cin>>D>>H>>B>>L;
        d=(B+0.0)/ceil((B+0.0)/D)/2.0;
        l=(L+0.0)/ceil((B+0.0)/D);
        double l1=0,r1=H;
        while(r1-l1>1e-5)
        {
            double m=(l1+r1)/2;
            if (Q(m)>0) r1=m;else l1=m;
        }
        if(cas>1) cout<<endl;
        printf("Case %d:
",cas);
        printf("%.2lf
",H-l1);

    }
    return 0;
}