UVALive-3485 Bridge(二分答案+自适应辛普森积分)

前言

本题是我在《训练指南》上发现的一道有趣的题目，有两种做法，一种是直接运用数学方法求积分，对这题而言较为繁琐；另一种则是运用自适应辛普森积分算法，简洁地求解此题。

《训练指南》一书中还有很多比较好的题目，以后我会抽空刷一刷，巩固自己的基础。

---

题目传送门

戳我，戳我

---

题目大意

你的任务是修建一座大桥。桥上等距地摆放着若干个塔，塔高为H，宽度忽略不计。相邻两座塔之间的距离不能超过D。塔之间的绳索形成全等的对称抛物线。桥长度为B，绳索总长为L，如下图所示求建最少的塔时绳索的最下端离地的高度y。

【输入格式】
输入第一行为测试数据组数T。每组数据包含4个整数D,H,B,L(B<=L)。

【输出格式】
对于每组数据，输出绳索底部离地高度，保留两位小数。

---

Solution

本题其实思考难度不高，但需要一些数学知识。
下面根据《训练指南》里的做法简要分析。

我们称相邻塔之间的部分为“间隔”。当塔数最小时，间隔数为$n=\lceil {B\over D}\rceil$，每个间隔宽度$D_1={B\over n}$，每段绳索长度为$L_1={L\over n}$。然后就可以根据$D_1$与$L_1$求出底部离地高度$y$。

由于$y$很难直接计算(显然)，所以考虑解方程。我们设宽度为$w$，高度为$h$的抛物线的长度为$p(w, h)$，我们只需要解一个关于$h$的方程$p(D_1,h)=L_1$，然后计算$y=H-h$。容易发现的是，$p(D_1,h)$关于$h$单调增，所以直接二分$h$。于是，问题只剩下$p(w,h)$该怎么求呢？

根据微积分知识，可导函数在区间$[a,b]$内的弧长为

$$\int_a^b\sqrt {1+[f'(x)]^2} dx$$
至于这个是怎么求出来的呢，有兴趣的可以看看下面简单的数学推算与说明，觉得显然或者无聊的跳过就行了。

我们在$x$轴上选取很小很小的一段，设其长度为$dx$，即$dx$->$0$。那么在$dx$上方的函数图像可以近似的看作一条线段，设线段在$y$轴上的投影长度为$dy$。那么，根据勾股定理，线段的长度就是$l=\sqrt{dx^2+dy^2}$。而${dy\over dx}=f'(x)$，所以我们将$dx$提出来，就有$l=\sqrt {1+[f'(x)]^2} dx$，我们求一整段区间的函数图像弧长就是$l$的积分。于是就有了上面的那条(莫名其)妙的公式。

由于每条抛物线都全等，所以函数解析式相同。我们将其中一个抛物线的顶点放在原点上，那么解析式就是$f(x)=ax^2$。我们又知道$w$和$h$，且$f(-{w \over 2}) = h$，解得$a={4h \over w^2}$。根据对称性，弧长为

$$2\int_0^{w\over2}\sqrt{1+[f'(x)]^2}dx=2\int_0^{w\over2}\sqrt{1+4a^2x^2}dx$$
求到这里，接下来我们只有两条路可以走了，第一条是求原函数的“不归路”，第二条是 辛普森积分的“康庄大道”。

我们先走”不归路”，试求不定积分。然而并不会。哦，根据定积分表里的公式，应该是这样子的

$$\int\sqrt{a^2+x^2}dx={1\over2}x\sqrt{a^2+x^2}+{1\over2}a^2ln|x+\sqrt{a^2+x^2}|+C$$
把上述公式变形为
$$2\int_0^{w\over2}\sqrt{1+4a^2x^2}dx=4a\int_0^{w\over2}\sqrt{({1\over {2a}})^2+x^2}dx$$
然后就可以套公式了。先贴代码：

---

Code(数学求积分)

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>
#define Eps 1e-5

using namespace std;

int T, D, H, B, L;

double F(double a, double x){
    double a2 = a * a, x2 = x * x;
    return (x * sqrt(a2 + x2) + a2 * log(fabs(x + sqrt(a2 + x2)))) / 2.0;
}

double Getlen(double w, double h){
    double a = 4.0 * h / (w * w);
    return 4.0 * a * (F(0.5/a, w/2.0) - F(0.5/a, 0.0));
}

int main(){

    freopen("LA3485.in", "r", stdin);
    freopen("LA3485.out", "w", stdout);

    scanf("%d", &T);

    for(int Case = 1; Case <= T; Case++){
      scanf("%d%d%d%d", &D, &H, &B, &L);

      int n = (B+D-1) / D;
      double D1 = (double)B / (double)n;
      double L1 = (double)L / (double)n;

      double low = 0.0, high = (double)H;
      while(low + Eps < high){
        double mid = (low + high) / 2.0;
        if(Getlen(D1, mid) - L1 < Eps)  low = mid;
        else  high = mid;
      }
      if(Case > 1)  printf("\n");
      printf("Case %d:\n%.2lf\n", Case, H-low);
    }


    return 0;
}

---

接下来就走simpson的道路。因为没有什么好写的，所以引用lrj的话，说明一下像辛普森积分这类方法的优点。

事实上，很多函数的积分都是没有封闭形式的，无法向上面那样推出公式以后再“套”。这种情况下，数值积分(numerical integral)是一个很好的工具。

---

Code(simpson)

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

int T, D, H, B, L;
double a;

double F(double x){
    return sqrt(1.0 + 4.0 * a * a * x * x);
}

double simpson(double L, double R){
    double mid = (L + R) / 2.0;
    return (F(L) + 4.0 * F(mid) + F(R)) * (R - L) / 6.0;
}

double integral(double L, double R, double Eps){
    double mid = (L + R) / 2.0;
    double ST = simpson(L, R), SL = simpson(L, mid), SR = simpson(mid, R);
    if(fabs(SL + SR - ST) <= 15.0 * Eps)  return SL + SR + (SL + SR - ST) / 15.0;
    return integral(L, mid, Eps/2.0) + integral(mid, R, Eps/2.0);
}


double Getlen(double w, double h){
    a = 4.0 * h / (w * w);
    return integral(0.0, w/2.0, 1e-5) * 2.0;
}

int main(){

    freopen("LA3485.in", "r", stdin);
    freopen("LA3485.out", "w", stdout);

    scanf("%d", &T);

    for(int Case = 1; Case <= T; Case++){
      scanf("%d%d%d%d", &D, &H, &B, &L);

      int n = (B+D-1) / D;
      double D1 = (double)B / (double)n;
      double L1 = (double)L / (double)n;

      double low = 0.0, high = (double)H;
      while(low + 1e-5 < high){
        double mid = (low + high) / 2.0;
        if(Getlen(D1, mid) - L1 < 1e-5)  low = mid;
        else  high = mid;
      }
      if(Case > 1)  printf("\n");
      printf("Case %d:\n%.2lf\n", Case, H-low);
    }


    return 0;
}

---

愿你有一天能和你重要的人重逢。