HDOJ-4498 Function Curve(自适应Simpson积分)

主要考点是数学计算和自适应Simpson积分

题意：

给你n系数，函数值，求在x=[0,100]这段区间内曲线的长度。

需要准备一些数学知识，高数/数分/工科数分/微积分 中曾经学过，曲线的长度=

可以确定，这是一个分段函数，想要计算积分值，就要确定每一段函数的表达式。要确定每一段函数的表达式，首先要确定每一段的转折点。

把y=100也算作一个函数，然后把抛物线的方程转为一下，变成： (kx^2) -(2akx) +(a^2k+b),x的二次，一次，零次系数为k,-2ak,a^2*k+b

这样得到n个抛物线（直线），再加上y=100这条直线，每两条线联立，解方程，算出交点，存到一个数组当中。最后把0和100也放到这里数组里当做起点和终点，排序。

然后遍历每一个小区间，扫描一遍所有系数得到使这个区间内函数最小的系数下标（考虑到y=100的情况在系数的数组里加入0,0,100），

然后对这个函数在小区间内进行自适应Simpson积分即可。

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
using namespace std;
const int N=100;
double A[N],B[N],C[N];
double X[N*N];
int tot,pos;

//被积函数，弧长公式
double F(double x)
{
    return sqrt(1.0+(x*2*A[pos]+B[pos])*(x*2*A[pos]+B[pos]));
}

double simpson(double a,double b)
{
    double c=a+(b-a)/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-a)/2;
    double L=simpson(a,c);
    double 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));
}

bool zero(double x)
{
    return (fabs(x)<1e-8);
}

/*解两条抛物线（直线）的交点
两方程联立(A1-A2)x^2+(B1-B2)x+(C1-C2)=0
相当于解方程 a*x^2+b*x+c=0*/
void cal(int a,int b,int c)
{
    double temp,diet;
    if(zero(a) && zero(b)) return;//无解
    if(zero(a))
    {
        temp=-c/b;
        if(temp>=0 && temp<=100)
            X[tot++]=temp;
        return;
    }
    diet=b*b-4*a*c;
    if(diet<0) return;        //无实数根
    else if(zero(diet))  //有一个实根
    {
        temp=-b/(2.0*a);
        if(temp>=0 && temp<=100)
            X[tot++]=temp;
        return;
    }
    else                      //有两个不同的实根
    {
        diet=fabs(diet);
        temp=(-b+sqrt(diet))/(2*a);
        if(temp>=0 && temp<=100)
            X[tot++]=temp;
        temp=(-b-sqrt(diet))/(2*a);
        if(temp>=0 && temp<=100)
            X[tot++]=temp;
    }
}

void work()
{
    int i,j,t,m,n;
    int a,b,k;
    double t1,t2,ans;

    memset(A,0,sizeof(A));
    memset(B,0,sizeof(B));
    memset(C,0,sizeof(C));
    memset(X,0,sizeof(X));
    ans=0;
    tot=0;
    A[0]=0;
    B[0]=0;
    C[0]=100;

    scanf("%d",&n);
    for(i=1;i<=n;i++)
    {
        //转为 A*x^2+b*x+k的形式,存储系数
        scanf("%d%d%d",&k,&a,&b);
        A[i]=k;
        B[i]=-2*a*k;
        C[i]=a*a*k+b;
    }

    //计算每两条抛物线的交点，并加入交点集合
    for(i=0;i<=n;i++)
        for(j=i+1;j<=n;j++)
            cal(A[i]-A[j],B[i]-B[j],C[i]-C[j]);
    X[tot++]=0;
    X[tot++]=100;
    sort(X,X+tot);
    //for(i=0;i<tot;i++) printf("%0.3lf\n",X[i]);

    //枚举每两个交点，并枚举找出两个交点之间的哪组系数使表达式最小
    for(i=0;i<tot-1;i++)
    {
        if(fabs(X[i]-X[i+1])<1e-8) continue;//区间为0
        double mid=(X[i]+X[i+1])/2.0;
        pos=0;//使表达式最小的系数下标
        t1=mid*mid*A[pos] + mid*B[pos] + C[pos];
        for(j=1;j<=n;j++)
        {
            t2=mid*mid*A[j] + mid*B[j] + C[j];
            if(t2<t1)
            {
                t1=t2;
                pos=j;
            }
        }
        ans+=asr(X[i],X[i+1],1e-8);
    }
    printf("%0.2lf\n",ans);

}

int main()
{
    //freopen("in.txt","r",stdin);
    int T;
    scanf("%d",&T);
    while(T--)
        work();
    return 0;
}