hdu4498 求曲线长度数值积分 simpson公式

最新推荐文章于 2022-12-27 17:35:35 发布

weixin_30329623

最新推荐文章于 2022-12-27 17:35:35 发布

阅读量294

点赞数

CC 4.0 BY-SA版权

文章标签： php

原文链接：http://www.cnblogs.com/--560/p/5329061.html

本文探讨了如何使用 Simpson 公式来估算连续曲线的定积分，特别关注了解二次方程的特殊情况。通过实例展示了在计算交点时的处理方式，以及如何利用二次方程的根来划分积分区间，进而应用 Simpson 方法进行精确计算。此方法在数值分析中具有重要应用价值。

http://acm.hdu.edu.cn/showproblem.php?pid=4498

simpson公式可以求连续曲线的定积分。

这里注意解ax^2+bx+c=0时a=0,退化为一次的情况。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<math.h>
#define REP(i,a,b) for(int i=a;i<=b;i++)
#define MS0(a) memset(a,0,sizeof(a))

using namespace std;

typedef long long ll;
const int maxn=1000100;
const int INF=1e9+10;
const double EPS=0.0000000001;

int n;
struct Line
{
    double k,a,b;
    void read()
    {
        scanf("%lf%lf%lf",&k,&a,&b);
    }
    double y(double x)
    {
        return k*(x-a)*(x-a)+b;
    }
};Line L[maxn];
struct Point
{
    double x,y;
    friend bool operator<(Point A,Point B)
    {
        return A.x<B.x;
    }
    void out()
    {
        printf("x=%.2f y=%.2f\n",x,y);
    }
};Point p[maxn];int pn;

bool can(double x,double y)
{
    if(x<=0||x>=100.0) return 0;
    REP(i,0,n){
        if(L[i].y(x)<y-EPS) return 0;
    }
    return 1;
}

void work(double a,double b,double c,Line A)
{
    double dt=b*b-4*a*c;
    if(dt<-EPS) return;
    if(fabs(a)<EPS){
        if(fabs(b)>EPS){
            double x=-c/b;
            if(can(x,A.y(x))) p[++pn]={x,A.y(x)};
        }
        return;
    }
    double x1=(-b-sqrt(dt))/(2*a);
    double x2=(-b+sqrt(dt))/(2*a);
    if(fabs(dt)<EPS){
        if(!can(x1,A.y(x1))) return;
        p[++pn]={x1,A.y(x1)};
    }
    else{
        if(can(x1,A.y(x1))) p[++pn]={x1,A.y(x1)};
        if(can(x2,A.y(x2))) p[++pn]={x2,A.y(x2)};
    }
}

void inter(Line A,Line B)
{
    double a=A.k-B.k;
    double b=-2*(A.k*A.a-B.k*B.a);
    double c=A.k*A.a*A.a-B.k*B.a*B.a+A.b-B.b;
    work(a,b,c,A);///解ax^2+bx+c=0，并取交点
}

void get_duan()
{
    double y=100.0;
    REP(i,1,n) y=min(y,L[i].y(0));
    p[++pn]={0,y};
    y=100.0;
    REP(i,1,n) y=min(y,L[i].y(100.0));
    p[++pn]={100.0,y};
}

void get_inter()
{
    pn=0;
    REP(i,0,n){
        REP(j,i+1,n){
            inter(L[i],L[j]);///取两条曲线的交点
        }
    }
    ///取端点0和100
    get_duan();
    sort(p+1,p+pn+1);
}

double F(double x,double k,double a,double b)
{
    double t=2*k*(x-a);
    return sqrt(1+t*t);
}
/// simpson求定积分
double simpson(double a,double b,double tk,double ta,double tb)
{
    double c=a+(b-a)/2;
    return (F(a,tk,ta,tb)+4*F(c,tk,ta,tb)+F(b,tk,ta,tb))*(b-a)/6;
}

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

double asr(double a,double b,double EPS,double tk,double ta,double tb)
{
    return asr(a,b,EPS,simpson(a,b,tk,ta,tb),tk,ta,tb);
}
///---
bool In(Line L,Point A)
{
    if(fabs(L.y(A.x)-A.y)<EPS) return 1;
    return 0;
}

double Len(Point A,Point B)
{
    if(fabs(A.x-B.x)<EPS) return 0;
    double mx=(A.x+B.x)/2;
    int cur=-1;
    REP(i,0,n){
        if(!In(L[i],A)||!In(L[i],B)) continue;
        if(cur==-1||L[i].y(mx)<L[cur].y(mx)) cur=i;
    }
    if(cur==-1) return 0;
    return asr(A.x,B.x,EPS,L[cur].k,L[cur].a,L[cur].b);
}

double solve()
{
    get_inter();///取交点
    /*
    REP(i,1,pn){
        p[i].out();
    }
    */
    double res=0;
    REP(i,1,pn-1){
        res+=Len(p[i],p[i+1]);///计算相邻两点的曲线长度
    }
    return res;
}

int main()
{
    freopen("in.txt","r",stdin);
    int T;cin>>T;
    while(T--){
        scanf("%d",&n);
        L[0].k=L[0].a=0;L[0].b=100.0;
        REP(i,1,n) L[i].read();
        printf("%.2f\n",solve());
    }
    return 0;
}