三分·三分求极值

最新推荐文章于 2021-01-17 18:57:46 发布

weixin_34092455

最新推荐文章于 2021-01-17 18:57:46 发布

阅读量280

点赞数

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/weierpeng/p/4393799.html

本文介绍使用三分法解决抛物线与点之间最短距离问题的方法。通过将问题转化为求凸函数极值，利用三分法逼近求解。文章详细解释了如何划分搜索区间，并给出了完整的C++代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

http://hihocoder.com/contest/hiho40/problem/1

描述

这一次我们就简单一点了，题目在此：

在直角坐标系中有一条抛物线y=ax^2+bx+c和一个点P(x,y)，求点P到抛物线的最短距离d。

提示：三分法

在之前的几周中我们了解到二分法作为分治中最常见的方法，适用于单调函数，逼近求解某点的值。
但当函数是凸形函数时，二分法就无法适用，这时就需要用到三分法。
从三分法的名字中我们可以猜到，三分法是对于需要逼近的区间做三等分：

我们发现lm这个点比rm要低，那么我们要找的最小点一定在[left,rm]之间。如果最低点在[rm,right]之间，就会出现在rm左右都有比他低的点，这显然是不可能的。同理，当rm比lm低时，最低点一定在[lm,right]的区间内。
利用这个性质，我们就可以在缩小区间的同时向目标点逼近，从而得到极值。

接下来我们回到题目上，抛物线和点之间的距离可以简单的用直线公式计算：
即d = min{sqrt((X - x)^2+(aX^2+bX+c-y)^2)}
该公式展开后为4次，需要采用求导等方法来求极值。对于计算机编程来说是很麻烦的一件事。

进一步观察题目，我们可以发现根据带入的X值不同，d的长度恰好满足凸形函数。
而我们要求的最短距离d，正好就是这个凸形函数的极值。
那么三分法不就正好可以用来解决这道题目了么？

需要注意在解题过程中一定要想清楚如何划分区间，我们求的各个变量到底是什么含义。
另外，这道题还有一个小小的trick，在解决的时候请一定要小心。

输入

第1行：5个整数a,b,c,x,y。前三个数构成抛物线的参数，后两个数x,y表示P点坐标。-200≤a,b,c,x,y≤200

输出

第1行：1个实数d，保留3位小数(四舍五入)

样例输入

2 8 2 -2 6

样例输出

2.437

 1 #include <iostream>
 2 #include <math.h>
 3 using namespace::std;
 4 int a;
 5 double Xp,Yp;
 6 const double precision=1E-6;
 7 
 8 inline double d(double x)
 9 {
10     return sqrt((Xp-x)*(Xp-x) + (Yp-a*x*x)*(Yp-a*x*x));
11 }
12 
13 
14 int main(int argc, const char* argv[])
15 {
16     int b,c;
17     int X0,Y0;
18     cin>>a>>b>>c>>X0>>Y0;
19     double D;
20     if (a==0) {
21         D=fabs(1.*(b*X0-Y0+c)/sqrt(b*b+1));
22     }else{
23         double Xc=-1.*b/(2*a);
24         double Yc=c-1.*b*b/(4*a);
25         Xp=X0-Xc;
26         Yp=Y0-Yc;
27         double x1,x2;
28         if (a<0){
29             a=-a;
30             Yp=-Yp;
31         }
32         if (Xp<0)Xp=-Xp;
33         if (Yp<0) {
34             x1=0.;
35             x2=Xp;
36         }else{
37             x1=Xp;
38             x2=sqrt(Yp/a);
39             if (x1>x2)swap(x1, x2);
40         }
41         double lmx,rmx,lm,rm;
42         while (fabs(d(x2)-d(x1))>precision) {
43             lmx=x1+1./3*(x2-x1);
44             lm=d(lmx);
45             rmx=x1+2./3*(x2-x1);
46             rm=d(rmx);
47             if (lm<rm){
48                 x2=rmx;
49             }else{
50                 x1=lmx;
51             }
52             
53         }
54         D=d((x1+x2)/2);
55     }
56     printf("%.3lf\n",D);
57     return 0;
58 }