LA 3485 (积分 辛普森自适应法) Bridge

桥的间隔数为n = ceil(B/D)，每段绳子的长度为L / n，相邻两塔之间的距离为 B / n

主要问题还是在于已知抛物线的开口宽度w 和 抛物线的高度h 求抛物线的长度

弧长积分公式为：

设抛物线方程为f(x) = ax²，则这段抛物线弧长为

查积分表或者自己分部积分算一下： 

二分抛物线高度x，使得每段抛物线长度为L / n，所求答案为H - x

#include <cstdio>
#include <cmath>

inline double F(double a, double x)
{//sqrt(a^2+x^2)的原函数
    double a2 = a*a, x2 = x*x;
    double s = sqrt(a2+x2);
    return (x*s + a2*log(x+s))/2;
}

double length(double w, double h)
{//宽为w，高为h的抛物线的长度
    double a = 4*h/w/w;
    double b = 0.5/a;
    return 4*a*(F(b, w/2) - F(b, 0));
}

int main()
{
    //freopen("in.txt", "r", stdin);

    int T;
    scanf("%d", &T);
    for(int kase = 1; kase <= T; kase++)
    {
        int D, H, B, L;
        scanf("%d%d%d%d", &D, &H, &B, &L);
        int n = (B-1)/D + 1;    //间隔数
        double d = (double)B / n; //间隔
        double l = (double)L / n;   //每段绳长
        double Left = 0, Right = H;
        while(Right - Left > 1e-5)
        {//二分求抛物线高度
            double mid = (Right + Left) / 2;
            if(length(d, mid) > l) Right = mid;
            else Left = mid;
        }
        if(kase > 1) puts("");
        printf("Case %d:
%.2f
", kase, H-Left);
    }

    return 0;
}

后面又介绍了一种Simpson自适应算法，可以求任意连续函数的积分。

维基百科辛普森法则

虽然不明白这个式子是怎么来的，但并不能阻止我们学习自适应辛普森算法。

书上还说可以将区间端点和中点的函数值作为参数传入以减少重复计算，求教。。

#include <cstdio>
#include <cmath>

double a;

inline double F(double x)
{ return sqrt(1+4*a*a*x*x); }

double simpson(double a, double b)
{
    double c = (a+b)/2;
    return (F(a)+4*F(c)+F(b))*(b-a)/6;
}

double asr(double a, double b, double eps, double A)
{
    double c = (a+b)/2;
    double L = simpson(a, c), R = simpson(c, b);
    if(fabs(L+R-A) <= 15*eps) return L + R + (L+R-A)/15.0;
    return asr(a, c, eps/2, L) + asr(c, b, eps/2, R);
}

double asr(double a, double b, double eps)
{
    return asr(a, b, eps, simpson(a, b));
}

double length(double w, double h)
{
    a = 4.0*h/w/w;
    return asr(0, w/2, 1e-5) * 2;
}

int main()
{
    //freopen("in.txt", "r", stdin);

    int T;
    scanf("%d", &T);
    for(int kase = 1; kase <= T; kase++)
    {
        int D, H, B, L;
        scanf("%d%d%d%d", &D, &H, &B, &L);
        int n = (B-1)/D + 1;
        double d = (double)B / n;
        double l = (double)L / n;
        double x = 0, y = H;
        while(y - x > 1e-5)
        {
            double m = (x + y) / 2;
            if(length(d, m) > l) y = m;
            else x = m;
        }
        if(kase > 1) puts("");
        printf("Case %d:
%.2f
", kase, H - x);
    }

    return 0;
}