Simpson积分学习笔记

昨天 ACM A C M 模拟的时候遇到了一道 Simpson S i m p s o n 积分相关的题

完全不知道怎么求,我们组 Fishmen F i s h m e n Bym B y m 嘲讽了很久

于是今天下午结合各种资料还是看了一下

这个东西不要觉得它看上去讲什么积分很高级 实际上认真推导也不是很难

首先 Simpson S i m p s o n 积分的式子是这个东西

这里写图片描述

积分号的下标代表区间左端点,上标代表区间右端点

这个值的意思就是函数在 [a,b] [ a , b ] 区间与 x x 轴,x=a x=b x = b 围成的面积

首先引入一些东西, Simpson S i m p s o n 积分就是用一个二次曲线来无限逼近原曲线

来达到求得近似值的效果

首先引入一个东西 对于一个二次函数 f(x)=ax2+bx+c f ( x ) = a x 2 + b x + c

求积分可得 F(x)=x0f(x)dx=a3x3+b2x2+cx+D F ( x ) = ∫ 0 x f ( x ) d x = a 3 x 3 + b 2 x 2 + c x + D

在这里 D D 是一个常数 这个东西具体怎么证等我问了物理组的同学再来填坑

那么LRf(x)dx=F(R)F(L)
RLf(x)dx=a3(R3L3)+b2(R2L2)+c(RL) ∫ L R f ( x ) d x = a 3 ( R 3 − L 3 ) + b 2 ( R 2 − L 2 ) + c ( R − L )
RLf(x)dx=(RL)[a3(L2+R2+LR)+b2(L+R)+c] ∫ L R f ( x ) d x = ( R − L ) [ a 3 ( L 2 + R 2 + L R ) + b 2 ( L + R ) + c ]
RLf(x)dx=ba6(2aL2+2aR2+2aLR+3bL+3bR+6c) ∫ L R f ( x ) d x = b − a 6 ( 2 a L 2 + 2 a R 2 + 2 a L R + 3 b L + 3 b R + 6 c )
RLf(x)dx=ba6{(aL2+bL+c)+(aR2+bR+c)+4[a(L+R2)2+b(L+R2)+c]} ∫ L R f ( x ) d x = b − a 6 { ( a L 2 + b L + c ) + ( a R 2 + b R + c ) + 4 [ a ( L + R 2 ) 2 + b ( L + R 2 ) + c ] }
RLf(x)dx=ba6[f(L)+f(R)+4f(L+R2)] ∫ L R f ( x ) d x = b − a 6 [ f ( L ) + f ( R ) + 4 f ( L + R 2 ) ]

就这样证完了… 然后我们来考虑这个代码怎么实现

double Simpson(double a, double b) {
    int mid = (a + b) / 2;
    return (b - a) / 6 * (f(a) + f(b) + 4 * f(mid));
}

好了写完了
这样子直接计算我们发现误差非常大,毕竟原图像可能不能很好的用二次函数逼近

自适应Simpson积分

那怎么办呢,我们可以考虑模拟微分的过程,把区间分成无穷小份再把得到的答案积分就是最后的答案了

但是一般我们只要求一个近似值,所以我们就递归到满足正常误差即可

还有一个问题 每次递归下去都有精度误差,那么累加起来可能就和正确答案相差甚远了

这个时候我们可以在递归的时候提高精度要求,尽量让小的答案精确这样就可以了

在刘汝佳的代码中判断是这样写的

double _asr(double l, double r, double v, double eps) {
    double mid = (l + r) / 2;
    double v1 = simpson(l, mid), v2 = simpson(mid, r);
    if(fabs(v1+v2-v) <= 15*eps) return v1+v2 + (v1+v2-v)/15;
    else return _asr(l, mid, v1, eps/2) + _asr(mid, r, v2, eps/2);
}

大概意思就是把多出的精度损失也加回去,这样子能够比较好的保证精确度

裸的模板题

这个就不用讲了 套板子即可

Codes
#include<bits/stdc++.h>

using namespace std;
typedef double db;

const db eps = 1e-6;
db a, b, c, d;

db f(db x) {
    return (c * x + d) / (a * x + b);
}

db integral(db x1, db x2) {
    db mid = (x1 + x2) / 2;
    return (x2 - x1) / 6 * (f(x1) + f(x2) + f(mid) * 4);
}

db Simpson(db x1, db x2, db eps) {
    db mid = (x1 + x2) / 2, res1, res2, res;
    res = integral(x1, x2); res1 = integral(x1, mid); res2 = integral(mid, x2);
    if(fabs(res1 + res2 - res) < eps * 15) return res1 + res2 + (res1 + res2 - res) / 15;   
    return Simpson(x1, mid, eps / 2) + Simpson(mid, x2, eps / 2);
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("4525.in", "r", stdin);
    freopen("4525.out", "w", stdout);
#endif
    db L, R;
    cin >> a >> b >> c >> d >> L >> R;
    printf("%.6lf", Simpson(L, R, eps));
    return 0;
}
  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值