数值分析：数值积分与数值微分

1 数值积分概述

1.1 引言

对于许多实际问题的求解往往需要计算积分。在高等数学中计算积分采用的是著名的牛顿--莱布尼兹公式： $\int_a^b{f(x)}dx=F(b)-F(a)$ 这里 $F(x)$ 是 $f(x)$ 的原函数。从理论上来说这个公式很完善，但是这个公式在实际应用中使用是很困难的。原因有3点：

求解原函数 $F(x)$ 往往是很困难的；
大部分积分即使求解出了原函数，其形式往往很复杂，计算过程中舍入误差很大；
当 $f(x)$ 是由测量或数值计算给出的数表时，牛顿--莱布尼兹公式也不能运用。

因此研究积分的数值计算方法是很有必要的。

1.2 数值积分基本方法

插值法
在积分区间 $[a,b]$ 上取一组点： $a\le x_0<x_1<...<x_n\le b$ ，构造函数 $f(x)$ 的 $n$ 次拉格朗日插值多项式为 $L_n(x)=\sum\limits_{k=0}^n{l_k(x)f(x_k)}$ 其中 $l_k(x)=\prod\limits_{i=0}^n{\frac{x-x_i}{x_k-x_i}}=\frac{\omega_{k+1}(x)}{(x-x_k)\omega'_n(x_k)}$ 这里 $\omega_{k+1}(x)=\prod\limits_{i=0}^n(x-x_i)$ 从而得到以下数值积分公式（机械求积公式）： $\int_a^bf(x)dx \approx \int_a^bL_n(x)dx=\sum\limits_{k=0}^nA_kf(x_k)$ 其中求积系数 $A_k$ 通过插值基函数 $l_k(x)$ 的积分得到，即 $A_k=\int_a^bl_k(x)dx$ 为插值型求积公式由第4章的拉格朗日插值余项定理可知，插值型的求积公式的余项为 $R[f]=\int_a^bf(x)dx-\int_a^bL_n(x)dx=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)dx$

1.3 代数精度（重要！！）

定义1 如果某个求积公式对于次数不超过 $m$ 的多项式均能准确地成立，但对于 $m+1$ 次多项式就不能准确成立，则称该求积公式具有 $m$ 次代数精度

求解方法：要使求积公式具有 $m$ 次代数精度，只要令它对于 $f(x)=1,x,...,x^m$ 能够准确成立，而对于 $f(x)=x^{m+1}$ 不能准确成立即可，公式如下：

$\left\{\begin{aligned}&\sum\limits_{k=0}^nA_k=b-a\\&\sum\limits_{k=0}^nA_kx_k=\frac{b^2-a^2}{2}\\& ... \\&\sum\limits_{k=0}^nA_kx_k^m=\frac{b^{m+1}-a^{m+1}}{m+1}\\\end{aligned}\right.$

定理1 对给定的 $n+1$ 个互异节点 $x_0,x_1,...,x_n\in [a,b]$ ，总存在求积系数 $A_k(k=0,1,...,n)$ 使得机械求积公式至少具有 $n$ 次代数精度

定理2 机械求积公式至少具有 $n$ 次代数精度的充要条件是它是插值型的。

1.4 求积公式的稳定性与收敛性

（暂略）

2 牛顿--柯特斯公式

2.1 牛顿-柯特斯公式的建立

求积节点在 $[a,b]$ 内等距分布式，插值型求积公式称为牛顿--柯特斯(Newton-Cotes)求积公式，下面给出具体形式：
在 $[a,b]$ 上取 $n+1$ 个等距节点 $x_k=a+kh(k=0,1,...,n)$ ，其中 $\frac{b-a}{n}$ ，令 $x=a+th$ 得到： $\int_a^bf(x)dx \approx (b-a)\sum\limits_{k=0}^nC_k^{(n)}f(x_k)$ 其中 $C_k^{(n)}=\frac{h}{b-a}\int_0^n\prod\limits_{j=0,j\neq k}^n\frac{t-j}{k-j}dt=\frac{(-1)^{n-k}}{nk!(n-k)!}\int_0^n\prod\limits_{j=0,j\neq k}^n(t-j)dt$ $C_k^{(n)}$ 称为柯特斯系数（可以通过查表获得）

柯特斯系数具有如下性质：

$\sum_{k=0}^nC_k^{(n)}=1$ ；
$C_{n-k}^{(n)}=C_k^{(n)}$

2.2 低阶牛顿-柯特斯公式

梯形公式
当 $n=1$ 时，牛顿-柯特斯公式变为 $\int_a^bf(x)dx \approx \frac{f(a)+f(b)}{2}(b-a)$
这就是梯形公式，容易验证其具有一次代数精度。
余项表示为：在 $[a,b]$ 内存在一点 $\eta$ ，使得 $R_T=\frac{f''(\eta)}{2!}\int_a^b(x-a)(x-b)dx=-\frac{f''(\eta)}{12}(b-a)^3$
辛普森(Simpson)公式
当 $n=2$ 时，牛顿-柯特斯公式变为 $\int_a^bf(x)dx \approx \frac{f(a)+4f(c)+f(b)}{6}(b-a)$ 其中 $c=\frac{a+b}{2}$ 这就是辛普森(Simpson)公式，容易验证它具有三次代数精度。
余项表示为：在 $[a,b]$ 内存在一点 $\eta$ ，使得 $R_S=\frac{f^{(4)}(\eta)}{4!}\int_a^b(x-a)(x-c)^2(x-b)dx=-\frac{f^{(4)}(\eta)}{2880}(b-a)^5$
柯特斯(Cotes)公式
当 $n=4$ 时，牛顿-柯特斯公式简称为柯特斯(Cotes)公式，具有五次代数精度。
余项表示为：在 $[a,b]$ 内存在一点 $\eta$ ，使得 $R_S=-2(b-a)\frac{f^{(6)}(\eta)}{924}(\frac{b-a}{4})^6$

2.3 复化牛顿-柯特斯公式

复化梯形公式
将 $[a,b]$ 区间 $n$ 等分，自区间长度 $h=\frac{b-a}{n}$ ，于是有复化梯形公式 $T_n=\frac{h}{2}\left[f(a)+2\sum\limits_{k=0}^{n-1}f(x_k)+f(b)\right]$ 其余项公式为:在 $[a,b]$ 内存在一点 $\eta$ $R_{T_n}=-\frac{h^2}{12}(b-a)f''(\eta)$
复化辛普森公式
设子区间 $[x_k,x_{k+1}]$ 的中点 $x_{k+\frac{1}{2}}$ ，复化辛普森公式为 $S_n=\frac{h}{6}\left[f(a)+4\sum\limits_{k=0}^{n-1}f(x_{k+\frac{1}{2}})+2\sum\limits_{k=1}^{n-1}f(x_k)+f(b)\right]$ 其余项公式为：在 $[a,b]$ 内存在一点 $\eta$ $R_{S_n}=-\frac{b-a}{180}(\frac{h}{2})^4f^{(4)}(\eta)$
算法设计（以复化辛普森公式为例）
（1）输入区间端点 $a,b$ 及扥分数 $\frac{n}{2}$ ( $n$ 为偶数)，半步长 $h=\frac{b-a}{n}$
（2）置 $f_i=f(a+ih)(i=1,2,...,n)$
（3）置 $P=f_1,Q=0$
（4）对 $j=2,3,...,n-2$ ，置 $P+f_{j+1}=>P,Q+f_j=>Q$
（5）输出 $S=\frac{h}{6}[f(a)+4P+2Q+f(b)]$

3 多重积分

考虑二重积分 $\iint\limits_R f(x,y)dA$ 是曲面 $z=f(x,y)$ 与平面区域 $R$ 围成的体积，对于矩形区域 $R=\{(x,y)|a \le x \le b,c\le y \le d \}$ 可写成累次积分 $\iint\limits_R f(x,y)dA=\int_a^b\left(\int_c^df(x,y)dy\right)dx$ 同样可以使用复化梯形公式与复化辛普森公式求解

复化梯形公式
将 $[a,b]$ ， $[b,c]$ 区间分为 $n$ ， $m$ 等份，步长 $h=\frac{b-a}{n}$ ， $k=\frac{d-c}{m}$ 先对积分 $\int_c^df(x,y)dy$ 应用复化梯形公式，则有 $\int_c^df(x,y)dy=\frac{k}{2}\left[f(x,c)+2\sum\limits_{i=0}^{n-1}f(x,y_i)+f(x,d)\right]$ 从而得到 $\int_a^b\left(\int_c^df(x,y)dy\right)dx=\frac{k}{2}\left[\int_a^bf(x,c)+2\sum\limits_{i=0}^{n-1}\int_a^bf(x,y_i)+\int_a^bf(x,d)\right]$ 再对每个积分分别应用复化梯形公式即可
复化辛普森公式
将 $[a,b]$ ， $[b,c]$ 区间分为 $n$ ， $m$ 等份，步长 $h=\frac{b-a}{n}$ ， $k=\frac{d-c}{m}$ 先对积分 $\int_c^df(x,y)dy$ 应用复化辛普森公式，令 $y_i=c+ik,y_{i+\frac{1}{2}}=c+(i+\frac{1}{2})$ 则有 $\int_c^df(x,y)dy=\frac{k}{6}\left[f(x,c)+4\sum\limits_{i=0}^{n-1}f(x,y_{i+\frac{1}{2}})+2\sum\limits_{i=1}^{n-1}f(x,y_i)+f(x,b)\right]$ 从而得到 $\int_a^b\left(\int_c^df(x,y)dy\right)dx=\frac{k}{6}\left[\int_a^bf(x,c)+4\sum\limits_{i=0}^{n-1}\int_a^bf(x,y_{i+\frac{1}{2}})+2\sum\limits_{i=1}^{n-1}\int_a^bf(x,y_i)+\int_a^bf(x,b)\right]$ 再对每个积分分别应用复化辛普森公式即可