数值积分知识
数值积分
对于增加插值节点序列: { x i } i = 0 n \left\{x_i\right\}_{i=0}^{n} {xi}i=0n,由插值定理给出:
f ( x ) = ∑ i = 0 n y i l i ( x ) + f ( n + 1 ) ( ξ ) ( n + 1 ) ! ∏ i = 0 n ( x − x i ) f(x)=\sum_{i=0}^{n}y_i l_i(x)+\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^{n}{(x-x_i)} f(x)=i=0∑nyili(x)+(n+1)!f(n+1)(ξ)i=0∏n(x−xi)
对 f ( x ) f(x) f(x)在区间 [ a , b ] [a,b] [a,b]上积分得到:
∫ a b f ( x ) = ∫ a b [ ∑ i = 0 n y i l i ( x ) + f ( n + 1 ) ( ξ ) ( n + 1 ) ! ∏ i = 0 n ( x − x i ) ] d x \int_{a}^{b}f(x)=\int_{a}^{b}\left[\sum_{i=0}^{n}y_i l_i(x)+\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^{n}{(x-x_i)}\right]\text{d}x ∫abf(x)=∫ab[i=0∑nyili(x)+(n+1)!f(n+1)(ξ)i=0∏n(x−xi)]dx
得到结论:
对于插值型的求积公式,其截断误差为:
∫ a b ( f ( n + 1 ) ( ξ ) ( n + 1 ) ! ∏ i = 0 n ( x − x i ) ) d x \boxed{\int_{a}^{b}\left(\frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^{n}{(x-x_i)}\right)\text{d}x} ∫ab((n+1)!f(n+1)(ξ)i=0∏n(x−xi))dx
其数值计算公式为:
∫ a b [ ∑ i = 0 n y i l i ( x ) ] d x \int_{a}^{b}\left[\sum_{i=0}^{n}y_i l_i(x)\right]\text{d}x ∫ab[i=0∑nyili(x)]dx
更特殊地
,我们对区间 [ a , b ] [a,b] [a,b]进行等间距插值选取步长 h = b − a n h=\dfrac{b-a}{n} h=nb−a,得到 n + 1 n+1 n+1个插值节点 x i = a + i h , i = 0 , 1 , 2 , 3 , . . . , n x_i=a+ih,i=0,1,2,3,...,n xi=a+ih,i=0,1,2,3,...,n,这就对著名的柯特斯公式
。
常用的数值求积公式
梯形公式
-
公式推导:
- 对于定积分 ∫ a b f ( x ) d x \int_{a}^{b}f(x)dx ∫abf(x)dx,将区间 [ a , b ] [a,b] [a,b]进行 n n n等分,每个小区间 [ x i , x i + 1 ] [x_{i},x_{i + 1}] [xi,xi+1]的长度 h = b − a n h=\frac{b - a}{n} h=nb−a,其中 x i = a + i h x_{i}=a+ih xi=a+ih, i = 0 , 1 , ⋯ , n i = 0,1,\cdots,n i=0,1,⋯,n。
- 在每个小区间 [ x i , x i + 1 ] [x_{i},x_{i + 1}] [xi,xi+1]上,用梯形面积近似小曲边梯形的面积。梯形面积公式为 = ( 上底 + 下底 ) × 高 2 =\frac{(上底 + 下底)\times高}{2} =2(上底+下底)×高,对于小区间 [ x i , x i + 1 ] [x_{i},x_{i + 1}] [xi,xi+1],上底为 f ( x i ) f(x_{i}) f(xi),下底为 f ( x i + 1 ) f(x_{i+1}) f(xi+1),高为 h h h。
- 则
∫ a b f ( x ) d x ≈ ∑ i = 0 n − 1 h 2 [ f ( x i ) + f ( x i + 1 ) ] = h 2 [ f ( x 0 ) + 2 ∑ i = 1 n − 1 f ( x i ) + f ( x n ) ] 。 \int_{a}^{b}f(x)dx\approx\sum_{i = 0}^{n - 1}\frac{h}{2}[f(x_{i})+f(x_{i + 1})]=\frac{h}{2}[f(x_{0}) + 2\sum_{i = 1}^{n - 1}f(x_{i})+f(x_{n})]。 ∫abf(x)dx≈i=0∑n−12h[f(xi)+f(xi+1)]=2h[f(x0)+2i=1∑n−1f(xi)+f(xn)]。
-
误差分析:
- 设 f ( x ) f(x) f(x)在 [ a , b ] [a,b] [a,b]上具有二阶连续导数,则梯形公式的误差 R T = − ( b − a ) 3 12 n 2 f ′ ′ ( ξ ) R_{T}=-\frac{(b - a)^{3}}{12n^{2}}f^{\prime\prime}(\xi) RT=−12n2(b−a)3f′′(ξ)
其中 ξ ∈ ( a , b ) \xi\in(a,b) ξ∈(a,b)。
- 设 f ( x ) f(x) f(x)在 [ a , b ] [a,b] [a,b]上具有二阶连续导数,则梯形公式的误差 R T = − ( b − a ) 3 12 n 2 f ′ ′ ( ξ ) R_{T}=-\frac{(b - a)^{3}}{12n^{2}}f^{\prime\prime}(\xi) RT=−12n2(b−a)3f′′(ξ)
辛普森公式
-
公式推导:
- 同样将区间 [ a , b ] [a,b] [a,b]进行 n n n等分( n n n为偶数),在每个子区间 [ x i , x i + 1 ] [x_{i},x_{i + 1}] [xi,xi+1]上,用二次函数 p ( x ) = A x 2 + B x + C p(x)=A x^{2}+Bx + C p(x)=Ax2+Bx+C来近似 f ( x ) f(x) f(x)。
- 辛普森公式为
∫ a b f ( x ) d x ≈ h 6 [ f ( x 0 ) + 4 ∑ i = 0 n − 1 f ( x k + 1 2 ) + 2 ∑ i = 1 n − 1 f ( x i ) + f ( x n ) ] \int_{a}^{b} f(x)\text{d}x\approx\frac{h}{6}\left[f(x_{0}) +4\sum\limits_{i = 0}^{n-1}f(x_{k+\frac{1}{2}})+2\sum_{i = 1}^{n - 1}f(x_{i})+f(x_{n})\right] ∫abf(x)dx≈6h[f(x0)+4i=0∑n−1f(xk+21)+2i=1∑n−1f(xi)+f(xn)]
其中 h = b − a n h = \frac{b - a}{n} h=nb−a
-
误差分析:
- 设 f ( x ) f(x) f(x)在 [ a , b ] [a,b] [a,b]上具有四阶连续导数,则辛普森公式的误差
R S = − ( b − a ) 5 180 ( 2 n ) 4 f ( 4 ) ( ξ ) R_{S}=-\frac{(b - a)^{5}}{180(2n)^{4}}f^{(4)}(\xi) RS=−180(2n)4(b−a)5f(4)(ξ)
- 设 f ( x ) f(x) f(x)在 [ a , b ] [a,b] [a,b]上具有四阶连续导数,则辛普森公式的误差
其中 ξ ∈ ( a , b ) \xi\in(a,b) ξ∈(a,b)
习题
例
对于函数 f ( x ) = sin x x f(x)=\frac{\sin x}{x} f(x)=xsinx,给出 n = 8 n = 8 n=8 时的函数表,试用复合梯形公式及复合辛普森公式计算积分
I = ∫ 0 1 sin x x d x I=\int_{0}^{1}\frac{\sin x}{x}dx I=∫01xsinxdx
并估计误差。
解
当 n = 8 n=8 n=8时,复合梯形公式得到:
T n = 0.125 2 [ f ( 0 ) + f ( 1 ) + 2 × ( ∑ i = 1 7 f ( x i ) ) ] = 0.9456909 T_n=\frac{0.125}{2}\left[f(0)+f(1)+2\times \left(\sum _{i=1}^{7}f(x_i)\right)\right]=0.9456909 Tn=20.125[f(0)+f(1)+2×(i=1∑7f(xi))]=0.9456909
其中 x i = i 8 x_i=\frac{i}{8} xi=8i
同理,复合辛普森公式得到:
S n = 0.25 6 [ f ( 0 ) + f ( 1 ) + 2 × ( f ( 0.25 ) + f ( 0.5 ) + f ( 0.75 ) ) + 4 × ( f ( 0.125 ) + f ( 0.325 ) + f ( 0.625 ) + f ( 0.875 ) ) ] = 0.9460832 S_n=\frac{0.25}{6}\left[f(0)+f(1)+2\times \left(f(0.25)+f(0.5)+f(0.75)\right)+4\times \left(f(0.125)+f(0.325)+f(0.625)+f(0.875)\right)\right]=0.9460832 Sn=60.25[f(0)+f(1)+2×(f(0.25)+f(0.5)+f(0.75))+4×(f(0.125)+f(0.325)+f(0.625)+f(0.875))]=0.9460832
为了求解 f ( x ) = sin x x f(x)=\frac{\sin x}{x} f(x)=xsinx的高阶导数:
sin x x = ∫ 0 1 cos ( x t ) d t \frac{\sin x}{x}=\int_{0}^{1}\cos(xt)\text{d}t xsinx=∫01cos(xt)dt
由维尔斯特拉斯定理不难证明,求导与积分可以换序:
f ( k ) ( x ) = ∫ 0 1 t k cos ( x t + k π 2 ) d x f^{(k)}(x)=\int_{0}^{1}t^{k}\cos \left( xt+\frac{k\pi}{2}\right)\text{d}x f(k)(x)=∫01tkcos(xt+2kπ)dx
其上界:
max f ( k ) ( x ) = max ∫ 0 1 t k cos ( x t + k π 2 ) d x ≤ 1 k + 1 \max f^{(k)}(x)=\max \int_{0}^{1}t^{k}\cos \left( xt+\frac{k\pi}{2}\right)\text{d}x\leq \frac{1}{k+1} maxf(k)(x)=max∫01tkcos(xt+2kπ)dx≤k+11
得到复合梯形公式误差:
∣ R 8 ( f ) ∣ = 1 12 ( 0.125 ) 2 × 1 3 \left|R_{8}(f)\right|=\frac{1}{12}\left(0.125\right)^2\times \frac{1}{3} ∣R8(f)∣=121(0.125)2×31
得到复合辛普森公式误差:
∣ R 4 ( f ) ∣ = 1 180 ( 0.25 2 ) 4 × 1 5 \left|R_{4}(f)\right|=\frac{1}{180}\left(\frac{0.25}{2}\right)^4\times \frac{1}{5} ∣R4(f)∣=1801(20.25)4×51