引子
有時候我們需要計算一個函式的定積分，粗略上可以使用估算的方法，如圖所示，將原本的曲線粗略地看成一個梯形，這個方法叫梯形法則（Trapezoidal Rule），也叫做一階牛頓-柯特斯閉型積分公式，
其中
所謂一階，指的就是n=1的情況，
最理想的情況就是把這個影像分割成 無數 個梯形，便可求出對應的定積分，
但是在實際操作的情況下，梯形法則為了保證速度無法取極多的點，這樣照成梯形法則誤差較大，
分割成無限個梯形其實就等效于
因此我們將考慮更高階的公式，本文將要介紹的便是二階牛頓-柯特斯閉型積分公式（辛普森法）， 即將函式近似看成一條拋物線，顯然一階牛頓-柯特斯閉型積分公式需要在首尾取兩個點方可得到f(x)的決議式，而二階要得到一個拋物線方程則需要取三個點，才能得到決議式，而牛頓-柯特斯閉型積分公式都采用等距取點的方法，所以辛普森法需要取首尾點以及中間點，
推導方法
通過辛普森法，只需要計算出邊界和中點的函式值就可以得到估算的函式的積分
綠色是被積函式，紅線是使用梯形法則得到的，藍線是辛普森法得到的，可以發現相對梯形法則，辛普森法相對而言更符合被積函式，繼續提
高n的值可以提高估算的準確度，但顯然也會提高計算的復雜性，而二階的對于圖中的函式基本上可以達到要求，
對于一段基本符合二次函式的函式，辛普森法可以得到較為準確的結果，然而當函式變得復雜時，辛普森將出現明顯的誤差，這時候該怎么計算積分，
擺在面前的有兩條路，一個是不斷細磁區間，再將區間求和，另一條是增加n的值，使得函式更加符合被積函式，本文主要討論第一種方法，
正文
走第一條路，如果采用將區間平均細分為若干份，求得各個區間的積分后再統一求和，這樣在計算量上和梯形法則上相等，但是精度提升不了多少，
然而，不得不承認，對于一個函式而言，存在二次擬合比較好的地方，也存在比較差的地方，如果在擬合比較好的地方細分就會浪費計算資源，而統一的平均細分，也導致在擬合比較差的地方達不到想要的精度，
因此，如何在擬合好的地方降低細分程度，擬合差的地方提高細分進度成了重中之重，
而本文想介紹的關鍵自適應辛普森積分法，就是解決這個問題的，
上式是自適應辛普森積分法是否對一個區間進行細分的判定條件，
其實是根據這個結論：
來源
原理其實很簡單，首先將這個區間進行二分為兩個區間肯定能得到更為精準的結果，如果這個結果和直接計算該區間的辛普森積分得到的結果差距大于某一個值，則認為對這個區間繼續細分是有意義的，反之差距小于該值則說明繼續細分對提高精度沒有什么幫助，就不進行下一步的細分，以降低計算的復雜度，這個值一般是esp的15倍（至于為什么是15倍可以看看這篇文章Notes on the Adaptive Simpson Quadrature Routine），調節esp可以控制估計值和精確值在多少位前是一樣的，
擴展：為什么二分一個區間能夠得到一個更精確的結果，從影像上很好理解，兩個函式要比一個函式擬合得更緊密，從數值上說明是因為辛普森演算法的誤差公式為
p是f(x)在[l,r]區間內使得f(x)四階導數取得最大值的值，具體證明可以參考A short proof of the error term in Simpson's rule（看不懂？本筆者也看不懂），根據琴生不等式可以推匯出區間一分為二后的誤差和小于不一分為二的誤差（大概可以證明吧，以后有時間證明一下QAQ），
據說（美）薩奧爾（Sauer, T.）著；裴玉茹，馬賡宇譯.數值分析（原書第2版），機械工業出版社，第240頁， 有相對完善的介紹，以后買一本看看，
代碼
#include <iostream> #include <iomanip> #include <cmath> using namespace std; double F(double x) { return x * log(x); } double PF(double x) { return 1.0 / 2.0 * x * x * log(x) - 1.0 / 4.0 * x * x; } //primitive function: 1/2xlnx-1/4x int count_simpson =0; double simpson(double l,double r) { count_simpson++; return (r - l) / 6.0 * (F(l) + F(r) + 4 * F((r + l) / 2)); } double asr( double l,//區間左邊界 double r,//區間右邊界 double integral,//該區間的辛普森積分 double esp ) { double mid = (l + r) / 2; double in_l = simpson(l, mid), in_r = simpson(mid, r); if (fabs(in_l + in_r - integral) > esp * 15) //細分有價值 return asr(l, mid, in_l, esp) + asr(mid, r, in_r, esp); else return in_l + in_r + (in_l + in_r - integral) / 15.0;//回傳這個肯定比回傳integral精度更高 } int main() { cout << setprecision(16)<<"估算" << asr(1, 8, simpson(1, 8), 1e-7) << " && 精確" << PF(8) - PF(1) << endl; cout << "計算次數:" << count_simpson << endl; return 0; }
運行結果：
擴展
為什么我非要用辛普森，這個演算法為什么不能用梯形法則代替辛普森的作用？其實是可以用梯形法則代替辛普森的，
但是無論是計算量還是精度來說效果都不是很好（QAQ），不用更高階的原因也很簡單，提高階增加編碼難度同時提高的精度難以彌補耗時，
后記
高考和強基計劃審核結束，筆者如愿以償地被理想的大學錄取，由于錄取專業與數學強相關，日后更新的檔案也將高度圍繞數學方面展開，這個排版調了我半天，emm，以后還是老老實實用markdown吧，
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/556707.html
標籤：其他
上一篇：云計算在商業運營中的潛力
下一篇：返回列表