辛普森积分

最新推荐文章于 2022-10-16 16:44:23 发布

psc233

最新推荐文章于 2022-10-16 16:44:23 发布

阅读量486

点赞数

分类专栏：辛普森文章标签：辛普森

本文链接：https://blog.csdn.net/qq_40875466/article/details/102458429

版权

辛普森专栏收录该内容

1 篇文章 0 订阅

订阅专栏

simpson积分公式

$\int_{a}^{b}f(x)\,dx\approx\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]$

他实质上就是用二次函数来拟合我们的函数
推导过程：
$g(x)=Ax^2+Bx+C$
$G(x)=\frac{Ax^3}{3}+\frac{Bx^2}{2}+Cx+D$
$\int_{a}^{b}f(x)\,dx=\int_{a}^{b}g(x)dx$
$= G (b) - G (a)$
$=\frac{Ab^3}{3}+\frac{Bb^2}{2}+Cb-\frac{Aa^3}{3}-\frac{Ba^2}{2}-Ca$
$=\frac{A}{3}(b-a)(a^2+ab+b^2)+\frac{B}{2}(b-a)(b+a)+C(b-a)$
$=\frac{b-a}{6}(2Aa^2+2Aab+2Ab^2+3Bb+3Ba+6C)$
$=\frac{b-a}{6}(Aa^2+Ba+C+Ab^2+Bb+C+(Aa^2+2Aab+Ab^2+2Bb+2Ba+4C))$
$=\frac{b-a}{6}(Aa^2+Ba+C+Ab^2+Bb+C+4(A(\frac{a+b}{2})^2+B(\frac{a+b}{2})+C))$
$=\frac{b-a}{6}(g(a)+g(b)+4g(\frac{a+b}{2}))$
我们在用simpson时用f函数去替代g函数即可
【模板】自适应辛普森法1

#include<bits/stdc++.h>
#define db double
using namespace std;
db a,b,c,d,l,r;
db f(db x){return (c*x+d)/(a*x+b);}
db calc(db l,db mid,db r){return (r-l)*(f(l)+4*f(mid)+f(r))/6;}
db simpson(db l,db r,db eps,db val){
	db mid=(l+r)/2.0;	
	db lmid=(l+mid)/2,rmid=(mid+r)/2;
	db vall=calc(l,lmid,mid),valr=calc(mid,rmid,r);
	if (fabs(vall+valr-val)<=eps) return val; else 
	return simpson(l,mid,eps/2.0,vall)+simpson(mid,r,eps/2.0,valr);
}
int main(){
	scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r);
	printf("%.6lf\n",simpson(l,r,1e-7,calc(l,(l+r)/2,r)));
}

圆的面积并

#pragma GCC optimize (3,"inline","Ofast")
#include<bits/stdc++.h>
#define db double
using namespace std;
const int MAXN=1010;
const db eps=1e-13;
struct cir{
    db x,y,r;
}o[MAXN],tmp[MAXN];
pair<db,db>po[MAXN];
int n,st,ed;
db xl[MAXN],xr[MAXN];
bool cmp(cir x,cir y)
{
    return x.r>y.r;
}
bool cmp2(cir x,cir y)
{
    return x.x-x.r<y.x-y.r;
}
double sqr(double x)
{
    return x*x;
}
double dis(cir x,cir y)
{
    sqr(x.x-y.x)+sqr(x.y-y.y);
}
bool pan(cir x,cir y)
{
    dis(x,y)<=sqr(x.r-y.r);
}
void cut()
{
    sort(o+1,o+n+1,cmp);
    int k=0,j;
    for (int i=1;i<=n;i++)
    {
        for (j=1;j<=k;j++)
        if (pan(o[i],tmp[j]))
         {
            break;
         }
        if (j>k) tmp[++k]=o[i];
    }
    n=k;
    for (int i=1;i<=n;i++) o[i]=tmp[i];
}
db get(db x)
{
    int tot=0;
    for (int i=st;i<=ed;i++)
    {
        if (xl[i]>=x||xr[i]<=x) continue;
        db dis=sqrt(o[i].r-sqr(x-o[i].x));
        po[++tot]=make_pair(o[i].y-dis,o[i].y+dis);
    }
    sort(po+1,po+tot+1);
    int j;
    db ans=0;
    for (int i=1;i<=tot;i++)
    {
        db l=po[i].first,r=po[i].second;
        for (j=i+1;j<=tot;j++)
         if (po[j].first<=r) r=max(r,po[j].second);
         else break;
        ans+=r-l;
        i=j-1;
    }
    return ans;
}
db calc(db x,db fl,db fmid,db fr) {return x/6.0*(fl+fmid*4+fr);}
db simpson(db l,db mid,db r,db fl,db fmid,db fr,db s)
{
    db mid1=(l+mid)/2,mid2=(mid+r)/2;
    db t1=get(mid1),t2=get(mid2);
    db g1=calc(mid-l,fl,t1,fmid),g2=calc(r-mid,fmid,t2,fr);
    if (fabs(g1+g2-s)<eps) return g1+g2;
    return simpson(l,mid1,mid,fl,t1,fmid,g1)+simpson(mid,mid2,r,fmid,t2,fr,g2);
}
void solve()
{
    cut();
    sort(o+1,o+n+1,cmp2);
    for (int i=1;i<=n;i++) xl[i]=o[i].x-o[i].r,xr[i]=o[i].x+o[i].r,o[i].r=sqr(o[i].r);
    int j;
    db ans=0;
    for (int i=1;i<=n;i++)
    {
        db l=xl[i],r=xr[i];
        for (j=i+1;j<=n;j++)
         if (xl[j]>r) break;
         else
         {
            r=max(r,xr[j]);
         }
        db mid=(l+r)/2;
        st=i; ed=j-1; i=j-1;
        db fl=get(l),fmid=get(mid),fr=get(r);
        ans+=simpson(l,mid,r,fl,fmid,fr,calc(r-l,fl,fmid,fr));
    }
    printf("%.3lf",ans);
}
int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++)
    {
        scanf("%lf%lf%lf",&o[i].x,&o[i].y,&o[i].r);
    }
    solve();
}

题目描述

小Q终于回家了。他开始了他最喜欢的Qazor之旅。Qazor的等离子场对过近的或过远的目标没有效果。小Q同时在多个位置释放了等离子场。他想知道他的等离子场能覆盖多大的面积。
输入
第一行一个整数n,表示等离子场的数量。
接下来n行,每行四个整数x,y,R,r,分别表示等离子场的横坐标,纵坐标,影响范围的上界，影响范围的下界。
输出
一行一个实数，保留三位小数，表示等离子场覆盖的面积。
样例输入
2
-1 0 2 1
1 0 2 1
样例输出
16.742
数据规模与约定
对于10%的数据,n=1
对于20%的数据,n<=2
对于40%的数据,n<=50
对于另外10%的数据,x=0,y=0
对于另外20%的数据,r=0
对于100%的数据,n<=200,|x|,|y|,|R|,|r|<=10000
时间限制 1s
内存限制 256M

#include<bits/stdc++.h>
#define db double
#define mk(x,y) make_pair(x,y)	
using namespace std;
const int inf=0x3f3f3f3f;
const int N=1010;
struct node{
	int x,y,R,r;
};
pair<db,int>c[N];
node v[N];
int n;
bool cmp(pair<db,int> x,pair<db,int> y){
	return x.first<y.first;
}
db sqr(db x){
	return x*x;
} 
db get(db x){
	int len=0;
	for (int i=1;i<=n;i++) {
	if (x>v[i].x-v[i].R&&x<v[i].x+v[i].R) {
		db d=sqrt(sqr(v[i].R)-sqr(v[i].x-x));
		c[++len]=mk(v[i].y-d,1); c[++len]=mk(v[i].y+d,-1);
	}
	if (x>v[i].x-v[i].r&&x<v[i].x+v[i].r) {
		db d=sqrt(sqr(v[i].r)-sqr(v[i].x-x));
		c[++len]=mk(v[i].y-d,-1); c[++len]=mk(v[i].y+d,1);
	} 
	}
	sort(c+1,c+len+1,cmp);
	int sum=1;
	db ans=0;
	for (int i=2;i<=len;i++) {
	if (sum) ans+=c[i].first-c[i-1].first;
	sum+=c[i].second;
	}
	return ans;
}
db calc(db l,db mid,db r){return (r-l)*(get(l)+4*get(mid)+get(r))/6;}
db simpson(db l,db r,db eps,db val){
	db mid=(l+r)/2;
	db lmid=(l+mid)/2,rmid=(mid+r)/2;
	db vall=calc(l,lmid,mid),valr=calc(mid,rmid,r);
	if (fabs(vall+valr-val)<eps) return val;
	else return simpson(l,mid,eps,vall)+simpson(mid,r,eps,valr); 
}
int main(){
	scanf("%d",&n);
	int lmin=inf,lmax=-inf;
	for (int i=1;i<=n;i++){
		scanf("%d%d%d%d",&v[i].x,&v[i].y,&v[i].R,&v[i].r);
		lmin=min(lmin,v[i].x-v[i].R); lmax=max(lmax,v[i].x+v[i].R); 
	}
	printf("%.3lf\n",simpson(lmin,lmax,1e-7,calc(1.0*lmin,1.0*(lmin+lmax)/2,1.0*lmax)));
}

psc233

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
辛普森积分

simpson积分公式∫abf(x) dx≈b−a6[f(a)+4f(a+b2)+f(b)]\int_{a}^{b}f(x)\,dx\approx\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]∫abf(x)dx≈6b−a[f(a)+4f(2a+b)+f(b)]他实质上就是用二次函数来拟合我们的函数推导过程：g(x)=Ax2+Bx+Cg(x)=...
复制链接

扫一扫