自適應Simpson法&&P4525 【模板】自適應辛普森法1
技術標籤：計算幾何
何為自適應Simpson法
自適應辛普森演算法(Adaptive Simpson′s rule)是一類近似演算法(Approximation algorithm)，主要用於在資訊計算時求較難反導的函式的定積分。其思想是利用二次函式曲線來不斷**擬合(Overfitting)所求曲線（顯然這傳統的直接用矩形或直邊梯形作為微元更精準），而所謂的Adapative（自適應）**則是用於優化時間複雜度的方法。 ----wikipedia
在計算幾何裡，顯然可在特殊情況下利用自適應Simpson法求圖形的面積。
Simpon公式及推導
Simpson公式
∫
a
b
f
(
x
)
d
x
≈
(
b
−
a
)
(
f
(
a
)
+
f
(
b
)
+
4
f
(
a
+
b
2
)
)
6
\int_a^bf(x)dx\approx\frac{(b-a)(f(a)+f(b)+4f(\frac{a+b}{2}))}{6}
推導
設 f ( x ) f(x) f(x)為被積函式， g ( x ) = A x 2 + B x + C g(x)=Ax^2+Bx+C g(x)=Ax2+Bx+C為用於擬合 f ( x ) f(x) f(x)的函式，即 g ( x ) ≈ f ( x ) g(x)\approx f(x) g(x)≈f(x)，則有：
∫
a
b
f
(
x
)
d
x
≈
∫
a
b
g
(
x
)
d
x
=
∫
a
b
(
A
x
2
+
B
x
+
C
)
d
x
=
A
3
(
b
3
−
a
3
)
+
B
2
(
b
2
−
a
2
)
+
C
(
b
−
a
)
\int_{a}^{b}f(x)dx\approx \int_{a}^{b}g(x)dx=\int_{a}^{b}(Ax^2+Bx+C)dx=\frac{A}{3}(b^3-a^3)+\frac{B}{2}(b^2-a^2)+C(b-a)
=
A
3
(
b
−
a
)
(
b
2
+
a
b
+
a
2
)
+
B
2
(
b
−
a
)
(
b
+
a
)
+
C
(
b
−
a
)
=
(
b
−
a
)
[
A
3
(
b
2
+
a
b
+
a
2
)
+
B
2
(
b
+
a
)
+
C
]
=\frac{A}{3}(b-a)(b^2+ab+a^2)+\frac{B}{2}(b-a)(b+a)+C(b-a)=(b-a)[\frac{A}{3}(b^2+ab+a^2)+\frac{B}{2}(b+a)+C]
= ( b − a ) 6 [ 2 A ( b 2 + a b + a 2 ) + 3 B ( b + a ) + 6 C ] =\frac{(b-a)}{6}[2A(b^2+ab+a^2)+3B(b+a)+6C] =6(b−a)[2A(b2+ab+a2)+3B(b+a)+6C]
= ( b − a ) 6 [ A a 2 + B a + C + A b 2 + B b + C + A a 2 + A b 2 + 2 A a b + 2 B a + 2 B b + 4 C ] =\frac{(b-a)}{6}[Aa^2+Ba+C+Ab^2+Bb+C+Aa^2+Ab^2+2Aab+2Ba+2Bb+4C] =6(b−a)[Aa2+Ba+C+Ab2+Bb+C+Aa2+Ab2+2Aab+2Ba+2Bb+4C]
= ( b − a ) 6 [ f ( a ) + f ( b ) + A ( a + b ) 2 + 2 B ( a + b ) + 4 C ] =\frac{(b-a)}{6}[f(a)+f(b)+A(a+b)^2+2B(a+b)+4C] =6(b−a)[f(a)+f(b)+A(a+b)2+2B(a+b)+4C]
= ( b − a ) 6 [ f ( a ) + f ( b ) + 4 A ( a + b 2 ) 2 + 4 B ( a + b 2 ) + 4 C ] =\frac{(b-a)}{6}[f(a)+f(b)+4A(\frac{a+b}{2})^2+4B(\frac{a+b}{2})+4C] =6(b−a)[f(a)+f(b)+4A(2a+b)2+4B(2a+b)+4C]
= ( b − a ) ( f ( a ) + f ( b ) + 4 f ( a + b 2 ) ) 6 =\frac{(b-a)(f(a)+f(b)+4f(\frac{a+b}{2}))}{6} =6(b−a)(f(a)+f(b)+4f(2a+b))
上述推導過程除了求了個一元二次多項式的積分外，運用的都是初等數學的知識。
除此推導方法之外，此公式還有一種更常規的推導方法，那就是利用牛頓-柯特斯求積公式。有篇相關文章寫的不錯，連結如下：
https://blog.csdn.net/qq_43069547/article/details/83153673
程式碼實現—自適應
如果想要利用Simpson公式求定積分，顯然我們要把積分割槽間分成多個足夠小的區間，但問題就在於如何去劃分才能既精確效率又高。
其實，Simpson法的程式實現很簡單，就是二分遞迴的過程，但關鍵在於“自適應”。所謂自適應，說的直白點，無非就是需要多分治幾次的地方，多分治幾次；不需要的則可以少分治幾次。
那麼，我們可以得出自適應Simpson法的核心實現過程：在二分遞迴劃分區間的過程中，如果滿足了精度需要，則可以終止當前遞迴，而這就是自適應辛普森法能夠自動控制區間分割大小的手段。
inline double queryans(double l,double r,double eps,double ans)//ans：利用Simpson法得到的以[l,r]為積分割槽間的f(x)的定積分； eps:精度
{
register double mid=(l+r)/2.0;
register double ansl=Simpson(l,mid),ansr=Simpson(mid,r);
if(fabs(ansl+ansr-ans)<=15.0*eps)//Problem 1
return ansl+ansr+(ansl+ansr-ans)/15.0;
return queryans(l,mid,eps/2.0,ansl)+queryans(mid,r,eps/2.0,ansr);
//Problem 2
}
上述程式碼有兩點需要說明的地方：
Problem 1：為什麼判斷是否還需要對當前區間繼續劃分擬合時的判斷條件是誤差 ≤ 15 \le15 ≤15倍的精度值？
答：據下面連結的文章說，wekipedia裡有相關證明，同時貼了張圖：
Problem 2：為什麼遞迴求折半積分割槽間的積分時精度？
如果要求當前區間 U U U計算的返回值誤差 e p s ( U ) < k eps(U)<k eps(U)<k的話，顯然其折半區間 L , R L,R L,R的返回值的誤差應有 e p s ( L ) ≤ k 2 , e p s ( R ) ≤ k 2 eps(L)\le \frac{k}{2},eps(R)\le \frac{k}{2} eps(L)≤2k,eps(R)≤2k這些不等關係，只有這樣才能保證 e p s ( U ) = e p s ( L ) + e p s ( R ) ≤ k 2 + k 2 = k eps(U)=eps(L)+eps(R)\le \frac{k}{2}+\frac{k}{2}=k eps(U)=eps(L)+eps(R)≤2k+2k=k。
最後，需要說明的是，本文大部分內容學習或借鑑自https://www.cnblogs.com/pks-t/p/10277958.html，這篇文章講的確實不錯，許多其他相關文章都沒提到的令人疑惑的地方點的較到位。
Code
#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;
double A,B,C,D,L,R;
inline double f(double x) { return (C*x+D)/(A*x+B); }
inline double Simpson(double a,double b) { return (f(a)+4*f((a+b)/2.0)+f(b))*(b-a)/6.0; }
inline double queryans(double l,double r,double eps,double ans)
{
register double mid=(l+r)/2.0;
register double ansl=Simpson(l,mid),ansr=Simpson(mid,r);
if(fabs(ansl+ansr-ans)<=15.0*eps) return ansl+ansr+(ansl+ansr-ans)/15.0;//15.0是科學推算得
return queryans(l,mid,eps/2.0,ansl)+queryans(mid,r,eps/2.0,ansr);//精度別忘除以2（列式等價得）
}
int main()
{
scanf("%lf%lf%lf%lf%lf%lf",&A,&B,&C,&D,&L,&R);
printf("%.6lf",queryans(L,R,1e-6,Simpson(L,R)));
return 0;
}