Simpson算法
用于求解定积分一类问题。
思路是用抛物线对原函数进行拟合,在误差足够小的前提下近似得到原函数的积分值。
Simpson算法是将原曲线近似成一段段抛物线后积分,取曲线上两点
A
(
a
,
0
)
,
B
(
b
,
0
)
A\left ( a, 0 \right ), B\left ( b, 0 \right )
A(a,0),B(b,0),然后等分,即
Δ
x
=
b
−
a
n
\Delta x = \frac{b - a}{n}
Δx=nb−a
对于每一段,作如下处理:
设区间起点为
x
2
i
−
1
x_{2i - 1}
x2i−1,终点为
x
2
i
x_{2i}
x2i,中点为
x
2
i
−
1
x_{2i - 1}
x2i−1,于是有利用过点
(
x
2
i
−
2
,
f
(
x
2
i
−
2
)
)
,
(
x
2
i
−
1
,
f
(
x
2
i
−
1
)
)
,
(
x
2
i
,
f
(
x
2
i
)
)
\left ( x_{2i - 2}, f\left ( x_{2i - 2} \right ) \right ), \left ( x_{2i - 1}, f\left ( x_{2i - 1} \right ) \right ), \left ( x_{2i}, f\left ( x_{2i} \right ) \right )
(x2i−2,f(x2i−2)),(x2i−1,f(x2i−1)),(x2i,f(x2i))的抛物线
g
(
x
)
=
A
x
2
+
B
x
+
C
g\left ( x \right ) = Ax^2 + Bx + C
g(x)=Ax2+Bx+C来取代
f
(
x
)
f\left ( x \right )
f(x),所以有:
f ( x 2 i − 2 ) = g ( x 2 i − 2 ) = A x 2 i − 2 2 + B x 2 i − 2 + C f\left ( x_{2i - 2} \right ) = g\left ( x_{2i - 2} \right ) = Ax_{2i - 2}^2 + Bx_{2i - 2} + C f(x2i−2)=g(x2i−2)=Ax2i−22+Bx2i−2+C
f ( x 2 i − 1 ) = g ( x 2 i − 1 ) = A x 2 i − 1 2 + B x 2 i − 1 + C = A ( x 2 i − 2 + x 2 i 2 ) 2 + B ( x 2 i − 2 + x 2 i 2 ) 2 + C f\left ( x_{2i - 1} \right ) = g\left ( x_{2i - 1} \right ) = Ax_{2i - 1}^2 + Bx_{2i - 1} + C = A\left ( \frac{x_{2i - 2} + x_{2i}}{2} \right )^2 + B\left ( \frac{x_{2i - 2} + x_{2i}}{2} \right )^2 + C f(x2i−1)=g(x2i−1)=Ax2i−12+Bx2i−1+C=A(2x2i−2+x2i)2+B(2x2i−2+x2i)2+C
f ( x 2 i ) = g ( x 2 i ) = A x 2 i 2 + B x 2 i + C f\left ( x_{2i} \right ) = g\left ( x_{2i} \right ) = Ax_{2i}^2 + Bx_{2i} + C f(x2i)=g(x2i)=Ax2i2+Bx2i+C
于是:
∫
x
2
i
−
2
x
2
i
f
(
x
)
d
x
≈
∫
x
2
i
−
2
x
2
i
g
(
x
)
d
x
=
(
A
2
x
3
+
B
2
x
2
+
C
x
)
∣
x
2
i
−
2
x
2
i
=
Δ
x
3
[
f
(
x
2
i
−
2
)
+
4
f
(
x
2
i
−
1
)
+
f
(
x
2
i
)
]
\begin{aligned} \int_{x_{2i - 2}}^{x_{2i}} f\left (x \right )\mathrm{d}x \approx & \int_{x_{2i - 2}}^{x_{2i}} g\left (x \right )\mathrm{d}x\\ = & \left ( \frac{A}{2} x^3 + \frac{B}{2}x^2 + Cx\right )\bigg|_{x_{2i - 2}}^{x_{2i}}\\ = & \frac{\Delta x}{3}\left [ f\left ( x_{2i - 2} \right ) + 4f\left ( x_{2i - 1} \right ) + f\left( x_{2i} \right ) \right ] \end{aligned}
∫x2i−2x2if(x)dx≈==∫x2i−2x2ig(x)dx(2Ax3+2Bx2+Cx)∣∣∣∣x2i−2x2i3Δx[f(x2i−2)+4f(x2i−1)+f(x2i)]
即
∫
a
b
f
(
x
)
d
x
≈
Δ
x
3
∑
i
=
0
2
n
−
2
[
f
(
x
2
i
)
+
4
f
(
x
2
i
+
1
)
+
f
(
x
2
i
+
2
)
]
\int_{a}^{b} f\left ( x \right ) \mathrm{d}x \approx \frac{\Delta x}{3} \sum_{i = 0}^{2n - 2} \left [ f \left( x_{2i} \right) +4f \left ( x_{2i + 1} \right ) + f \left ( x_{2i + 2} \right ) \right ]
∫abf(x)dx≈3Δxi=0∑2n−2[f(x2i)+4f(x2i+1)+f(x2i+2)]
特别地,当
n
=
1
n = 1
n=1时,上式转化为
∫
a
b
f
(
x
)
d
x
≈
b
−
a
6
[
f
(
a
)
+
4
f
(
a
+
b
2
)
+
f
(
b
)
]
\int_{a}{b} f \left ( x \right ) \mathrm{d}x \approx \frac{b - a}{6} \left [ f \left ( a \right ) + 4f \left ( \frac{a + b}{2} \right ) + f \left ( b \right ) \right ]
∫abf(x)dx≈6b−a[f(a)+4f(2a+b)+f(b)]
称为三点Simpson法。
自适应Simpson算法
自适应Simpson算法在三点Simpson算法的基础上提高了精度和效率:
- 对区间 [ a , b ] \left [ a, b \right ] [a,b]取中点 m i d = a + b − a 2 mid = a + \frac{b - a}{2} mid=a+2b−a;
- 分别对区间 [ a , b ] \left [ a, b \right ] [a,b]、区间 [ a , m i d ] \left [ a, mid \right ] [a,mid]和区间 [ m i d , b ] \left [ mid, b \right ] [mid,b]运用三点Simpson算法;
- 当 S 0 S_0 S0与 S 1 + S 2 S_1 + S_2 S1+S2的误差小于 e p s eps eps时,可近似得到曲线积分,如若不然,则继续递归调用。
模板如下:
inline double function(double t)
{
//所求积分函数
}
inline double simpson(double m, double n)
{
double mid = m + (n - m) / 2;
return (n - m) / 6 * (function(m) + 4 * function(mid) + function(n)); //拟合后进行数学整理得到的面积大小
}
double asr(double a, double b, double eps, double s)
{
double mid = a + (b - a) / 2;
double ls = simpson(a, mid), rs = simpson(mid, b);
if (fabs(ls + rs - s) < 16 * eps) //中点左右两部分的面积和与整个面积相差很小则可得到积分值
return ls + rs + (ls + rs - s) / 16;
return asr(a, mid, eps / 2, ls) + asr(mid, b, eps / 2, rs);
}