參考

• [學習筆記]自適應辛普森（Simpson）積分 [ 1 ] [1] [1]

引入

【模板】自適應辛普森法2 | 洛谷 P4526

• 給定 a a a 為小于等于 50 50 50 的實數，求
  ∫ 0 ∞ x a x ? x d x \int_0^\infin x^{\frac{a}{x}-x}dx ∫0∞?xxa??xdx
  若發散，輸出 o r z orz orz，
  若收斂，保留 5 5 5 位小數，

? \lceil ?辛普森法求積分 ? \rfloor ?

原理

• 對于 ∫ a b f ( x ) d x \int_a^bf(x)dx ∫ab?f(x)dx，如果這個 f ( x ) f(x) f(x) 我們很難對其進行積分計算的話
  我們可以用二次函式擬合，即設 f ( x ) ≈ ( A x 2 + B x + C ) f(x)\approx(Ax^2+Bx+C) f(x)≈(Ax2+Bx+C)
  然后經過推導 [ 1 ] [1] [1] ，得到一下結果：
  ∫ a b f ( x ) d x ≈ ( b ? a ) ( f ( a ) + 4 f ( a + b 2 ) + f ( b ) ) 6 \int_a^b f(x)dx\approx\cfrac{(b-a)(f(a)+4f(\frac{a+b}{2})+f(b))}{6} ∫ab?f(x)dx≈6(b?a)(f(a)+4f(2a+b?)+f(b))?
  其中 b ? a b-a b?a 的值越小，值越精確，

思考

• 但是，原函式 f ( x ) f(x) f(x) 不一定能用一個二次函式來得到很好的擬合，
  我們想，能否 f ( x ) f(x) f(x) 分成一段一段的，每一段都可以很好地用二次函式來擬合？
  貌似可行，但是怎么分段？這就引出了自適應辛普森法，

? \lceil ?自適應辛普森法求積分 ? \rfloor ?

原理

• 一個很自然的想法是，利用分治去處理區間，
  （1）分治 a s r ( l , r , e p s , a n s ) asr(l,r,eps,ans) asr(l,r,eps,ans)
  l , r l,r l,r 表示這個分治的區間的坐標
  e p s eps eps 表示一個期望令人滿意的精度
  a n s ans ans 存的是利用辛普森法算出的該區域的近似結果
  （2）進行分治計算，
  算出中點坐標 m i d = l + r 2 mid=\frac{l+r}{2} mid=2l+r?
  用辛普森法計算
  左半區域的近似和 l s u m = s i m p ( l , m i d ) lsum=simp(l,mid) lsum=simp(l,mid) 和右半區域的近似和 r s u m = s i m p ( m i d , r ) rsum=simp(mid,r) rsum=simp(mid,r)
  如果該區域的近似和與左半區域的近似和加上右半區域的近似和差不多，那么我們回傳 a n s ans ans，
  否則，繼續左遞回和右遞回，
• 一些奇怪的優化
  根據歷來 d a l a o dalao dalao 的優化，該區域精度需要有一個奇妙的判斷，但這不影響我們的理解，
  為什么遞回條件為精度*15

優點

• 可以計算某個函式的定積分，不需要考慮到如何去積分，
  該函式只要平滑，都可以帶入去運算，
  一些難以運算但可以通過近似去解的函式

缺點

• 首先，因為是實數運算，精度 和傳入 e p s eps eps 掛鉤，如果精度高了，運算時間 難免會增加，
  其次，要求函式較平滑，否則二次函式擬合程度較差，
  還有，要求積分上下限不能太大，更不能積到無窮去，時間和精度都不允許，
• 總結：精度差，時間高，函式要求平滑，上下限不能差太大，

代碼

double func(double x){
	/// 這里是 f(x) 的內容
    return sin(x);
}
double simp(double l,double r){		/// simpson 積分
    double mid = (l+r)/2.0;
    return (func(l) + 4.0*func(mid) + func(r)) * (r - l) / 6.0;
}
double asr(double l,double r,double eps,double ans){		/// 分治計算
    double mid = (l+r)/2.0;
    double ls = simp(l,mid),rs = simp(mid,r);
    if(fabs(ls + rs - ans) < 15 * eps)return ls + rs + (ls + rs - ans) / 15;
    return asr(l,mid,eps/2.0,ls) + asr(mid,r,eps/2.0,rs);
}

例題

模板：洛谷 P4526

∫ 0 ∞ x a x ? x d x \int_0^\infin x^{\frac{a}{x}-x}dx ∫0∞?xxa??xdx

• 什么時候發散呢？ a < 0 a<0 a<0發散，（打表也可以知道）
• 積分上下限怎么選擇？
  （1）積分下限不能取 0 0 0 ，因為該點計算無意義，下限取 e p s eps eps
  （2）積分上限：正無窮？做不了呀？畫圖：
  
  可以看到，最猛的 a = 50 a=50 a=50，也在 x > 10 x>10 x>10 出很快趨于 0 0 0了，
  我們積分上限取 15 15 15 即可，

double a;
double func(double x){
    return pow(x,a/x-x);
}
double simp(double l,double r){
    double mid = (l+r)/2;
    return (func(l) + 4*func(mid) + func(r)) * (r - l) / 6;
}
double asr(double l,double r,double eps,double ans){
    double mid = (l+r)/2;
    double ls = simp(l,mid),rs = simp(mid,r);
    if(fabs(ls + rs - ans) < 15 * eps)return ls + rs + (ls + rs - ans) / 15;
    return asr(l,mid,eps/2,ls) + asr(mid,r,eps/2,rs);
}
int main()
{
    cin >> a ;
    if(a < 0){
        puts("orz");
        return 0;
    }
    cout << fixed << setprecision(5) << asr(EPS,15,EPS,0);
    return 0;
}

橢圓：HDU1724

• Ellipse | HDU1724
  給定橢圓曲線 x 2 a 2 + y 2 b 2 = 1 \cfrac{x^2}{a^2}+\cfrac{y^2}{b^2}=1 a2x2?+b2y2?=1
  給定兩根直線 x = L x=L x=L 和 x = R x=R x=R
  求兩條直線之內，和橢圓所圍出的面積，
• 算出 y = f ( x ) y=f(x) y=f(x) 即可， f ( x ) = b 1 ? x 2 a 2 f(x)=b\sqrt{1-\cfrac{x^2}{a^2}} f(x)=b1?a2x2? ?
  然后帶入即可，

double a,b,l,r;
double func(double x){
    return sqrt(1.0-x*x/a/a)*b;
}
double simp(double l,double r){
    double mid = (l+r)/2.0;
    return (func(l) + 4.0*func(mid) + func(r)) * (r - l) / 6.0;
}
double asr(double l,double r,double eps,double ans){
    double mid = (l+r)/2.0;
    double ls = simp(l,mid),rs = simp(mid,r);
    if(fabs(ls + rs - ans) < 15 * eps)return ls + rs + (ls + rs - ans) / 15;
    return asr(l,mid,eps/2.0,ls) + asr(mid,r,eps/2.0,rs);
}
int main()
{
    int T;scanf("%d",&T);
    while(T--){
        scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
        printf("%.3f\n",asr(l,r,EPS,0)*2);
    }
    return 0;
}

Youmu with Master spark：QLU contest 003

• 二月八號新鮮的比賽，也是這題讓我了解到了自適應辛普森積分，
  但是這題出題人事先不知道辛普森積分做法 ???
• 【題意】Youmu with Master spark
  給定一個積分函式
  f ( x ) = ∫ 0 1 t 2 ∣ e 2 t ? x ∣ d t f(x)=\int_0^1t^2|e^{2t}-x|dt f(x)=∫01?t2∣e2t?x∣dt
  求
  ∑ i = 1 n f ( i ) \sum_{i=1}^n f(i) i=1∑n?f(i)
• 【資料范圍】
  樣例組數 T ≤ 1 0 5 T\le 10^5 T≤105
  1 ≤ n ≤ 1 0 9 1\le n\le 10^9 1≤n≤109
  精度誤差要求 ≤ 1 0 ? 8 \le 10^{-8} ≤10?8 （但是根據評測機顯示不需要精度這么高？）
• 【思路】
  看一下那個絕對值，當 x ≥ e 2 t x\ge e^{2t} x≥e2t 恒成立的話貌似可以直接去掉了，下限最大值為 e 2 e^2 e2
  也就是說，如果 x ≤ 7 x\le 7 x≤7，函式比較復雜，我們用自適應辛普森法去算出每一個 f ( i ) f(i) f(i) 來，
  如果 x ≥ 8 x\ge 8 x≥8，我們可以把函式進行化簡：
  f ( x ≥ 8 ) = ∫ 0 1 t 2 ( x ? e 2 t ) d t = x ∫ 0 1 t 2 d t ? ∫ 0 1 t 2 e 2 t d t = x 3 ? 1 2 ∫ 0 1 t 2 e 2 t d ( 2 t ) = x 3 ? 1 2 ∫ 0 2 ( y 2 ) 2 e y d y = x 3 ? 1 8 ∫ 0 2 y 2 d ( e y ) = x 3 ? 1 8 ( y 2 e y ∣ 0 2 ? ∫ 0 2 e y d ( y 2 ) ) = x 3 ? 1 8 ( 4 e 2 ? 2 ∫ 0 2 y d ( e y ) ) = x 3 ? 1 8 ( 4 e 2 ? 2 ( y e y ∣ 0 2 ? ∫ 0 2 e y d y ) ) = x 3 ? 1 8 ( 4 e 2 ? 2 ( 2 e 2 ? e y ∣ 0 2 ) ) = x 3 ? 1 8 ( 4 e 2 ? 2 ( 2 e 2 ? e 2 + 1 ) ) = x 3 ? 1 8 ( 2 e 2 ? 2 ) = x 3 ? e 2 ? 1 4 \begin{aligned} f(x\ge8)&=\int_0^1t^2(x-e^{2t})dt\\ &=x\int_0^1t^2dt-\int_0^1t^2e^{2t}dt\\ &=\frac{x}{3}-\frac{1}{2}\int_0^1t^2e^{2t}d(2t)\\ &=\frac{x}{3}-\frac{1}{2}\int_0^2(\frac{y}{2})^2e^ydy\\ &=\frac{x}{3}-\frac{1}{8}\int_0^2y^2d(e^y)\\ &=\frac{x}{3}-\frac{1}{8}\Big(y^2e^y\Big|_0^2-\int_0^2 e^yd(y^2)\Big)\\ &=\frac{x}{3}-\frac{1}{8}\Big(4e^2-2\int_0^2yd(e^y)\Big)\\ &=\frac{x}{3}-\frac{1}{8}\Big(4e^2-2(ye^y\Big|_0^2-\int_0^2e^ydy)\Big)\\ &=\frac{x}{3}-\frac{1}{8}\Big(4e^2-2(2e^2-e^y\Big|_0^2)\Big)\\ &=\frac{x}{3}-\frac{1}{8}\Big(4e^2-2(2e^2-e^2+1)\Big)\\ &=\frac{x}{3}-\frac{1}{8}(2e^2-2)\\ &=\frac{x}{3}-\frac{e^2-1}{4} \end{aligned} f(x≥8)?=∫01?t2(x?e2t)dt=x∫01?t2dt?∫01?t2e2tdt=3x??21?∫01?t2e2td(2t)=3x??21?∫02?(2y?)2eydy=3x??81?∫02?y2d(ey)=3x??81?(y2ey∣∣∣?02??∫02?eyd(y2))=3x??81?(4e2?2∫02?yd(ey))=3x??81?(4e2?2(yey∣∣∣?02??∫02?eydy))=3x??81?(4e2?2(2e2?ey∣∣∣?02?))=3x??81?(4e2?2(2e2?e2+1))=3x??81?(2e2?2)=3x??4e2?1??
  好家伙！主要是要用一次變數替換和兩次分部積分，復習了一下，
  那么， ? x ≥ 8 \forall x\ge 8 ?x≥8， f ( x ) = x 3 ? e 2 ? 1 4 f(x)=\frac{x}{3}-\frac{e^2-1}{4} f(x)=3x??4e2?1?，我們簡單做一下合并：
  A n s ( n ≥ 8 ) = ∑ i = 1 7 f ( i ) + ∑ i = 8 n f ( i ) = ∑ i = 1 7 f ( i ) + ∑ i = 8 n ( x 3 ? e 2 ? 1 4 ) Ans(n\ge8)=\sum_{i=1}^7f(i)+\sum_{i=8}^nf(i)=\sum_{i=1}^7f(i)+\sum_{i=8}^n\Big(\frac{x}{3}-\frac{e^2-1}{4}\Big) Ans(n≥8)=i=1∑7?f(i)+i=8∑n?f(i)=i=1∑7?f(i)+i=8∑n?(3x??4e2?1?)
  右邊的東西明顯可以化簡嘛！化簡就是最終答案了：
  A n s ( n ≥ 8 ) = ∑ i = 1 7 f ( i ) + ( 8 + n ) ( n ? 7 ) 2 × 3 + e 2 ? 1 4 ( n ? 7 ) Ans(n\ge 8)=\sum_{i=1}^7f(i)+\cfrac{(8+n)(n-7)}{2\times3}+\frac{e^2-1}{4}(n-7) Ans(n≥8)=i=1∑7?f(i)+2×3(8+n)(n?7)?+4e2?1?(n?7)

/*
 _            __   __          _          _
| |           \ \ / /         | |        (_)
| |__  _   _   \ V /__ _ _ __ | |     ___ _
| '_ \| | | |   \ // _` | '_ \| |    / _ \ |
| |_) | |_| |   | | (_| | | | | |___|  __/ |
|_.__/ \__, |   \_/\__,_|_| |_\_____/\___|_|
        __/ |
       |___/
*/
double ans[10];
double func(double x){
    return x*x*fabs(exp(2.0*x) - a);
}
double simp(double l,double r){
    double mid = (l+r)/2.0;
    return (func(l) + 4.0*func(mid) + func(r)) * (r - l) / 6.0;
}
double asr(double l,double r,double eps,double ans){
    double mid = (l+r)/2.0;
    double ls = simp(l,mid),rs = simp(mid,r);
    if(fabs(ls + rs - ans) < 15 * eps)return ls + rs + (ls + rs - ans) / 15;
    return asr(l,mid,eps/2.0,ls) + asr(mid,r,eps/2.0,rs);
}
int main()
{
    for(int i = 1;i <= 7;++i){
        a = i;
        ans[i] = asr(0,1,EPS,0);
    }
    const double ex = (exp(2.0) - 1)/4.0;
    int T;scanf("%d",&T);
    while(T--){
        int n;scanf("%d",&n);
        double res = 0;
        for(int i = 1;i <= min(n,7);++i)
            res += ans[i];
        if(n > 7){
            res += 1.0*(8.0+n)*(n-7.0)/6.0;
            res -= ex * (n - 7);
        }
        printf("%.10f\n",res);
    }
    return 0;
}