積分之自適應辛普森法 [學習筆記]

~hsm~發表於2019-03-25
筆記

文章目錄

宣告

本文證明和積分知識普及基本來自xyz32768。

一.預備技能

積分

積分是微積分學與數學分析裡的一個核心概念。通常分為定積分和不定積分兩種。直觀地說，對於一個給定的正實值函式，在一個實數區間上的定積分可以理解為在座標平面上，由曲線、直線以及軸圍成的曲邊梯形的面積值（一種確定的實數值）。
積分的一個嚴格的數學定義由波恩哈德·黎曼給出（參見條目“黎曼積分”）。黎曼的定義運用了極限的概念，把曲邊梯形設想為一系列矩形組合的極限。從十九世紀起，更高階的積分定義逐漸出現，有了對各種積分域上的各種型別的函式的積分。比如說，路徑積分是多元函式的積分，積分的區間不再是一條線段（區間[a,b]），而是一條平面上或空間中的曲線段；在面積積分中，曲線被三維空間中的一個曲面代替。對微分形式的積分是微分幾何中的基本概念。——摘自某度

其實，沒那麼麻煩，這樣解釋，可能更容易懂些（就是加顏色的話）：
積分（integral）的幾何意義是函式的曲線上 xx 的一段區間與 xx 軸圍成的曲邊梯形的面積：
在這裡插入圖片描述
xx 的區間為[a,b][a,b]
公式：∫abf(x)dx\int ^{b}_{a}f\left( x\right) dx

計算方法一：分割成無窮多個小區間。
∫abf(x)dx=lim⁡x→∞∑i=1nb−anf(a+b−ani)\int ^{b}_{a}f\left( x\right) dx=\lim _{x\rightarrow \infty }\sum ^{n}_{i=1}\dfrac {b-a}{n}f\left(a +\dfrac {b-a}{n}i\right)

計算方法二：牛頓-萊布尼茨公式。
定義:
如果函式f(x)f(x)在區間[a,b][a,b]上連續，並且存在原函式F(x)F(x) ，
則
∫abf(x)dx=F(b)−F(a)=F(x)∣ab\int ^{b}_{a}f\left( x\right) dx=F(b)-F(a)=F\left( x\right) | ^{b}_{a}

如果容易求出 nn 趨近於無窮大時 ff 的和，可以使用方法一，如 f(x)=x2f(x)=x^{2}。
而這時f(x)=1xf(x)=\dfrac {1}{x}就不適用。
如果容易求得 FF ，可以使用方法二，如 f(x)=1xf(x)=\dfrac {1}{x} 。
但如果兩個特點都不滿足，那麼兩種方法都無法使用。
於是，我們引入了數值積分。
最常用的就是自適應辛普森積分。

二.辛普森公式

借用一下大佬xyz32768的證明：
∫abf(x)dx≈∫ab(Ax2+Bx+C)dx\int ^{b}_{a}f\left( x\right) dx\approx\int ^{b}_{a}\left( Ax^{2}+Bx+C\right) dx
=A3(b3−a3)+B2(b2−a2)+C(a−b)=\frac{A}{3}(b^3-a^3)+\frac{B}{2}(b^2-a^2)+C(a-b)
=(b−a)6(2A(b2+ab+a2)+3B(b+a)+6C)=\frac{(b-a)}{6}(2A(b^2+ab+a^2)+3B(b+a)+6C)
=(b−a)6(2Ab2+2Aab+2Aa2+3Bb+3Ba+6C)=\frac{(b-a)}{6}(2Ab^2+2Aab+2Aa^2+3Bb+3Ba+6C)
=(b−a)6(Aa2+Ba+C+Ab2+Bb+C+4A(a+b2)2+4B(a+b2)+4C)=\frac{(b-a)}{6}(Aa^2+Ba+C+Ab^2+Bb+C+4A(\frac{a+b}{2})^2+4B(\frac{a+b}{2})+4C)
=(b−a)6(f(a)+f(b)+4f(a+b2))=\frac{(b-a)}{6}(f(a)+f(b)+4f(\frac{a+b}{2}))
我們得到辛普森公式：
∫abf(x)dx≈(b−a)(f(a)+f(b)+4f(a+b2))6\int ^{b}_{a}f\left( x\right) dx\approx \dfrac {\left( b-a\right) \left( f\left( a\right) +f\left( b\right) +4f\left( \dfrac {a+b}{2}\right) \right) }{6}
b−a 的值越小，上式兩邊越接近。

三.處理精度問題

同樣來自上面的大佬。

有了 SimpsonSimpson 公式，一個自然的想法是把積分割槽間拆成多個小區間後求和。
但是分成區間的個數和長度因積分割槽間和精度要求甚至被積函式而異。
換句話說，分的區間數太少不滿足精度要求，太多了會 TLE 。
「自適應」就是求積分時能夠自動控制切割的區間大小和長度，使精度能滿足要求。
具體地： solve(l,r,f)solve(l,r,f)表示求 ∫lrf(x)dx\int ^{r}_{l}f(x)dx。
（1）取 mid=l+r2mid=\dfrac {l+r}{2}。
（2）用 SimpsonSimpson 公式近似計算 f(x)f(x) 在區間 [l,mid][l,mid]和區間 [mid,r][mid,r]內的積分，分別為 lsls 和 rsrs ，及 [l,r][l,r]的近似積分 sumsum 。
（3）如果 ls+rsls+rs 與 sumsum的誤差允許，則返回 sumsum 。
（4）否則遞迴 solve(l,mid)solve(l,mid)和 solve(mid,r)solve(mid,r)，返回這兩個遞迴計算結果的和。

四.例題

P4525 【模板】自適應辛普森法1

title

LUOGU 4525

題目描述
計算積分
∫LRcx+dax+bdx\int ^{R}_{L}\dfrac {cx+d}{ax+b}dx
結果保留至小數點後6位。
資料保證計算過程中分母不為0且積分能夠收斂。
輸入輸出格式
輸入格式：
一行，包含6個實數a,b,c,d,L,R
輸出格式：
一行，積分值，保留至小數點後6位。
輸入輸出樣例
輸入樣例#1：
1 2 3 4 5 6
輸出樣例#1：
2.732937
說明
a,b,c,d∈[-10,10]
-100≤L<R≤100 且 R-L≥1

analysis

嗯，感覺是在套板子寫。

code

#include<bits/stdc++.h>
const double eps=1e-12;
static double a,b,c,d;

inline double F(double x)
{
    return (c*x+d)/(a*x+b);
}

inline double Simpson(double l,double r)
{
    return (r-l)*(F(l)+4.0*F((l+r)/2.0)+F(r))/6.0;
}

inline double Adaptive(double l,double r,double ans)
{
    double mid=(l+r)/2.0,left=Simpson(l,mid),right=Simpson(mid,r);
    if (fabs(left+right-ans)<eps)
        return left+right;
    else
        return Adaptive(l,mid,left)+Adaptive(mid,r,right);
}

int main()
{
    static double L,R;
    scanf("%lf %lf %lf %lf %lf %lf",&a,&b,&c,&d,&L,&R);
    printf("%.6lf",Adaptive(L,R,Simpson(L,R)));
    return 0;
}

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

title

LUOGU 4526

題目描述
計算積分
∫0∞xax−xdx\int ^{\infty }_{0}x^{\dfrac {a}{x}-x}dx
保留至小數點後5位。若積分發散，請輸出"orz"。
輸入輸出格式
輸入格式：
一行，包含一個實數，為a的值
輸出格式：
一行，積分值或orz
輸入輸出樣例
輸入樣例#1：
2.33
輸出樣例#1：
1.51068
說明
a<=50
請注意時空限制。

analysis

這道同樣是套板子，不過不要被式子的無窮大嚇到了。將函式影象畫出來，可以輕易地發現，當大於一定值時，函式趨於00，對答案基本沒有貢獻。這個值隨著aa增大增大，然而當a=50a=50時，這個值不超過88，為了求穩就可以從epseps積分到2020。為什麼不從00開始?因為當a=0a=0時，f(0)=1f(0)=1而不是00，所以計算會出現錯誤。避免特殊點對積分產生影響也是辛普森自適應積分的一個重點。記得epseps開小一點，最好開小44個數量級。

什麼時候無解?畫一下圖可以發現，當a&lt;0a&lt;0時，lim⁡x→0f(x)→+∞\lim\limits_{x\to 0} f(x)\to+\infty，函式不收斂。直接判掉即可。——Great_Influence

code

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-9;
static double a;

inline double F(double x)
{
    return pow(x,(a/x)-x);
}//原函式

inline double Simpson(double l,double r)
{
    return (r-l)*(F(l)+4.0*F((l+r)/2.0)+F(r))/6.0;
}//三點辛普森法

inline double Adaptive(double l,double r,double ans)
{
    double mid=(l+r)/2.0,left=Simpson(l,mid),right=Simpson(mid,r);
    if (fabs(left+right-ans)<=eps)
    	return ans;
    else
    	return Adaptive(l,mid,left)+Adaptive(mid,r,right);
}//自適應辛普森法

int main()
{
    scanf("%lf",&a);
    if (a<0)
        return puts("orz"),0;//判斷收斂/發散
    printf("%.5lf",Adaptive(eps,20,Simpson(eps,20)));
    return 0;
}

summary

程式碼真心不長，證明略有難度，不過記住結論就行了，然後，很重要的是，告訴我們，積分並不難（除了不定積分，這玩意兒我到現在都不知道是啥，只知道有個這東西）。

相關文章