为了提高数值积分的精确度,常将区间 [ a , b ] [a,b] [a,b] 等分成 n n n 个子区间
- 其长度为 h = b − a n h=\frac {b-a}n h=nb−a,
- 分点为 x k = a + k h ( k = 0 , 1 , . . . , n ) x_k = a + kh(k=0,1,...,n) xk=a+kh(k=0,1,...,n),
- 在每个子区间上用低阶的求积公式,然后将所有子区间上的计算结果加起来,这样得出的公式称为复化求积公式
这里介绍两个求积公式,复化梯形公式和复化辛普森公式
复化梯形公式
在每个子区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk,xk+1],用梯形公式
∫ a b f ( x ) d x ≈ b − a 2 [ f ( a ) + f ( b ) ] \int_{a}^b f(x)dx\approx \frac {b-a}2[f(a)+f(b)] ∫abf(x)dx≈2b−a[f(a)+f(b)]
得
T n = 1 n ∑ k = 0 n − 1 b − a 2 [ f ( x k ) + f ( x k + 1 ) ] T_n = \frac 1n \sum_{k=0}^{n-1} \frac {b-a}2 [ f(x_k) + f(x_{k+1}) ] Tn=n1k=0∑n−12b−a[f(xk)+f(xk+1)]
由
h
=
b
−
a
n
x
0
=
a
,
x
n
=
b
∑
k
=
0
n
−
1
f
(
x
k
)
=
f
(
x
0
)
+
f
(
x
1
)
+
⋯
+
f
(
x
n
−
1
)
∑
k
=
0
n
−
1
f
(
x
k
+
1
)
=
f
(
x
1
)
+
f
(
x
2
)
+
⋯
+
f
(
x
n
)
h=\frac {b-a}n \\ x_0=a,x_n=b\\ \sum_{k=0}^{n-1}f(x_k) = f(x_0)+f(x_1)+\cdots+f(x_{n-1}) \\ \sum_{k=0}^{n-1}f(x_{k+1}) = f(x_1)+f(x_2)+\cdots+f(x_n)
h=nb−ax0=a,xn=bk=0∑n−1f(xk)=f(x0)+f(x1)+⋯+f(xn−1)k=0∑n−1f(xk+1)=f(x1)+f(x2)+⋯+f(xn)
展开化简得
T n = h 2 [ f ( a ) + 2 ∑ k = 1 n − 1 f ( x k ) + f ( b ) ] T_n = \frac h2[ f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b) ] Tn=2h[f(a)+2k=1∑n−1f(xk)+f(b)]
依据上式,使用循环求和,即可实现复化梯形公式
代码如下:
// 复化梯形公式
double CompositeTrapezoidal(double a,double b,double n,double(*f)(double x))
{
double h = (b - a) / n;
double ans = f(a) + f(b);
for(int i = 1;i < n;i++)
{
ans += 2*f(a+i*h); //xk = a + kh ,
}
ans *= h/2;
return ans;
}
复化辛普森公式
记字段
[
x
k
,
x
k
+
1
]
[x_k,x_{k+1}]
[xk,xk+1]的中点为
x
k
+
1
2
=
x
k
+
x
k
+
1
2
=
a
+
(
2
k
+
1
)
h
2
x_{k+\frac 12}=\frac {x_k+x_{k+1}}2 =a+\frac {(2k+1)h}2
xk+21=2xk+xk+1=a+2(2k+1)h ,
在每个子区间
[
x
k
,
x
k
+
1
]
[x_k,x_{k+1}]
[xk,xk+1],用辛普森(Simpson)公式
∫ a b f ( x ) d x ≈ b − a 6 [ f ( a ) + 4 f ( a + b 2 ) + f ( b ) ] \int_{a}^b f(x)dx\approx \frac {b-a}6[f(a)+4f(\frac {a+b}2)+f(b)] ∫abf(x)dx≈6b−a[f(a)+4f(2a+b)+f(b)]
得
S n = ∑ k = 0 n − 1 h 6 [ f ( x k ) + 4 f ( x k + 1 2 ) + f ( x k + 1 ) ] S_n = \sum_{k=0}^{n-1} \frac h6 [f(x_k)+4f(x_{k+ \frac 12})+f(x_{k+1})] Sn=k=0∑n−16h[f(xk)+4f(xk+21)+f(xk+1)]
按上文的方式展开化简得
S n = h 6 [ f ( a ) + 4 ∑ k = 0 n − 1 f ( x k + 1 2 ) + 2 ∑ k = 1 n − 1 f ( x k ) + f ( b ) ] S_n=\frac h6[f(a)+4\sum_{k=0}^{n-1} f(x_{k+ \frac 12})+2\sum_{k=1}^{n-1}f(x_k)+f(b)] Sn=6h[f(a)+4k=0∑n−1f(xk+21)+2k=1∑n−1f(xk)+f(b)]
令
∑
k
=
0
n
−
1
f
(
x
k
+
1
2
)
=
∑
k
=
1
n
f
(
x
k
−
1
2
)
2
∑
k
=
1
n
−
1
f
(
x
k
)
=
2
∑
k
=
1
n
f
(
x
k
)
−
2
f
(
x
n
)
=
2
∑
k
=
1
n
f
(
x
k
)
−
2
f
(
b
)
\sum_{k=0}^{n-1} f(x_{k+ \frac 12})=\sum_{k=1}^{n} f(x_{k- \frac 12})\\ 2\sum_{k=1}^{n-1}f(x_k)=2\sum_{k=1}^{n}f(x_k)-2f(x_n)=2\sum_{k=1}^{n}f(x_k)-2f(b)\\
k=0∑n−1f(xk+21)=k=1∑nf(xk−21)2k=1∑n−1f(xk)=2k=1∑nf(xk)−2f(xn)=2k=1∑nf(xk)−2f(b)
最后整理得
S n = h 6 { f ( a ) − f ( b ) + ∑ k = 1 n [ 4 f ( x k − 1 2 ) + 2 f ( x k ) ] } S_n=\frac h6 \{f(a)-f(b)+\sum_{k=1}^{n}[ 4f(x_{k- \frac 12})+2f(x_k)] \} Sn=6h{f(a)−f(b)+k=1∑n[4f(xk−21)+2f(xk)]}
易知 x k − 1 2 = a + ( 2 k − 1 ) h 2 x_{k-\frac 12} =a+\frac {(2k-1)h}2 xk−21=a+2(2k−1)h,复化辛普森公式的实现与复化梯形公式类似
代码如下:
double CompositeSimpson(double a,double b,double n,double(*f)(double x))
{
double h = (b - a) / n;
double ans = f(a) - f(b);
double x = a;
for(int i=1;i<=n;i++)
{
ans += 4*f(a+(2*i-1)/2.0*h); // x(k-1/2)
ans += 2*f(a+i*h); // xk
}
ans *= h/6;
return ans;
}
测试
- 用复化梯形公式、复化辛普森公式求积分
∫
1
2
2
3
x
3
e
x
2
d
x
=
e
4
\int_1^2 \frac 23 x^3 e^{x^2}dx=e^4
∫1232x3ex2dx=e4,
取n=10
完整代码
#include <iostream>
#include <cmath>
using namespace std;
// 复化梯形公式
double CompositeTrapezoidal(double a,double b,double n,double(*f)(double x))
{
double h = (b - a) / n;
double ans = f(a) + f(b);
for(int i = 1;i < n;i++)
{
ans += 2*f(a+i*h);
}
ans *= h/2;
return ans;
}
// 复化辛普森公式
double CompositeSimpson(double a,double b,double n,double(*f)(double x))
{
double h = (b - a) / n;
double ans = f(a) - f(b);
double x = a;
for(int i=1;i<=n;i++)
{
ans += 4*f(a+(2*i-1)/2.0*h);
ans += 2*f(a+i*h);
}
ans *= h/6;
return ans;
}
int main()
{
double a,b,h;
int n;
cout<<"输入a :";cin>>a;
cout<<"输入b :";cin>>b;
cout<<"输入n :";cin>>n;
auto f=[](double x){return (2.0/3*x*x*x*exp(x*x));};
cout<<"精确值"<< exp(4)<< endl;
cout<<"使用复化梯形公式求得的结果:"<<CompositeTrapezoidal(a,b,n,f)
<<",绝对误差为: "<<abs(CompositeTrapezoidal(a,b,n,f)-exp(4))<<endl;
cout<<"使用复化辛普森公式求得的结果: "<<CompositeSimpson(a,b,n,f)
<<",绝对误差为: "<<abs(CompositeSimpson(a,b,n,f)-exp(4))<<endl;
return 0;
}
运行结果