Simpson 公式
[int_l^rf(x)dxapprox frac{r-l}{6}[f(l)+4f(frac{l+r}{2})+f(r)]
]
Simpson 公式将 ((l,f(l)),(r,f(r)),(frac{l+r}{2},f(frac{l+r}{2}))) 看作一个抛物线,求得近似解
具体过程如下:
令 (g(x)=ax^2+bx+c) 过 ((l,f(l)),(r,f(r)),(frac{l+r}{2},f(frac{l+r}{2}))) 三点
积分求得
[G(x)=int_0^xg(x)dx= frac{1}{3}ax^3+frac{1}{2}bx^2+cx+d
]
那么有
[egin{eqnarray*}
&& int_l^rg(x)dx=G(r)-G(l)\
&& = frac{1}{3}a(r^3-l^3)+frac{1}{2}b(r^2-l^2)+c(r-l)\
&& = frac{r-l}{6}[2a(l^2+lcdot r+r^2)+3b(l+r)+6c]\
&& = frac{r-l}{6}{(acdot l^2+bcdot l+c)+(acdot r^2+bcdot r+c)+4[a(frac{l+r}{2})^2+b(frac{l+r}{2})+c]}\
&& = frac{r-l}{6}[g(l)+4g(frac{l+r}{2})+g(r)]
end{eqnarray*}]
对于一个抛物线能取到等号,而将 (f(x)) 在区间 ([l,r]) 的端点和中点近似看作一个抛物线,就能取到一个近似值
double f(double x) { return ...; }
double simpson(double x) { return (r-l)*(f(l)+4*f((l+r)/2)+f(r))/6; }
自适应辛普森法
求出对区间 ([l,r]) 使用 Simpson 公式的结果,再求出对区间 ([l,mid]) 和 ([mid,r]) 使用的结果之和,若两者误差在精度要求范围内,就直接返回结果,否则递归求解两个子区间
double asr(double l,double r,double ans) {
double mid = (l+r)/2,L = simpson(l,mid),R = simpson(mid,r);
if (fabs(L+R-ans) <= eps) return L+R;
return asr(l,mid,L)+asr(mid,r,R);
}