当前位置：首页 > news >正文

数值分析—数值积分

news 来源：原创 2025/7/13 7:24:12

研究背景

积分的数学解法为牛顿莱布尼兹公式，数学表示为 $\int_{a}^{b} f(x)dx=F(b)-F(a)$ ，但应用该方法有如下困难：
1， $f (x)$ 的原函数有时不能用初等函数表示，如 $e^{-x^2}$ 和 $\frac{sinx}{x}$ ；
2，原函数可以用初等函数表示，但很复杂，如 $x^2\sqrt{1+2x^2}$ ；
3， $f (x)$ 本身无表达式，只有在若干点上的值。

针对这些问题，本章介绍几种近似求解积分的方法。

数值积分公式

数值积分

由中值定理 $\int_{a}^{b} f(x)dx=f(\xi)(b-a)，$ 其中 $f(\xi)$ 为 $f (x)$ 在[a,b]上的平均高度，可以不求原函数计算积分，但 $\xi$ 未知，故基于该方法，提出数值积分公式——用a代替 $\xi$ ，即 $\int_{a}^{b} f(x)dx\approx f(a)(b-a)$ ，计算左端点为宽的矩阵公式，称为左矩阵公式，同理还有右矩阵公式 $f (b) (b - a)$ 和中矩阵公式 $f(\frac{a+b}{2})(b-a)$ 。

矩形的误差可能过于大了，该方法则再次修正为 $\frac{f(a)+f(b)}{2}(b-a)$ ，即根据端点函数值计算梯形面积近似积分值。

该方法的再次改进为近似计算函数的平均值， $f(\xi)(b-a)\approx \frac{\sum_{k=0}^{n}f(x_{k})}{n+1}(b-a)$ ，用多个点的平均值代替 $f(\xi)$ ，上式可抽象为 $\sum_{k=0}^{n}A_{k}f(x_{k})$ ，其中 $A_{k}$ 称为权，该算法关键在于如何确定 $A_{k}$ 和 $x_{k}$ ，取点越多一定越精确，但如何衡量精确程度呢？

误差分析 $R(x)=\int_{a}^{b} f(x)dx-\sum_{k=0}^{n}A_{k}f(x_{k})$

代数精度

定义：若数值积分公式对一切 $\le m$ 次的多项式均能精确成立，则称次求积公式至少具有m次代数精度，若求积公式对一切 $\le m$ 次的多项式均能精确成立，但对 $m + 1$ 次多项式不能精确成立，则称此公式的代数精度为m次。

求积公式代数精度越高，使公式精确成立的多项式次数越高。

定理：
求积公式具有m次代数精度的充要条件，当 $f(x)=1,x,x^{2}...x^{m}$ 时公式均能准确成立，但 $f(x)=x^{m+1}$ 时不能准确成立。

例题：梯形公式的代数精度为多少？
解：梯形公式 $\int_{a}^{b} f(x)dx=\frac{f(a)+f(b)}{2}(b-a)$
令 $f (x) = 1$ ，左边=b-a=右边
令 $f (x) = x$ ，左= $\frac{b^{2}-a^{2}}{2}$ =右
令 $f(x)=x^{2}$ ， $\frac{b^3-a^3}{2}!=\frac{(a^2+b^2)}{2}(b-a)$ =右边
故该公式的代数精度为1。

插值型求积公式及Newton-Cotes公式

插值型求积公式

利用Lagrange插值求n次多项式的 $L_{n}(x)=\sum_{k=0}^{n}f(x_{k})l_{k}(x)$ ，
其中 $f(x)=L_{n}(x)+R_{n}(x),R_n{x}=\frac{f^{n+1}(\xi)}{(n+1)!}w_{n+1}(x)$
则 $\int_{a}^{b} f(x)dx=\int_{a}^{b} L_{n}(x)dx+\int_{a}^{b}R_{n}(x)dx=\int_{a}^{b} \sum_{k=0}^{n}f(x_{k})l_{k}(x)dx+\int_{a}^{b}R_{n}(x)dx=\sum_{k=0}^{n}f(x_{k})\frac{\int_{a}^{b}l_{k}(x)dx+\int_{a}^{b}R_{n}(x)dx}{A_{k}}$
即用插值多项式代替被积函数。

若定义求积公式 $\int_{a}^{b} f(x)dx\approx\sum_{k=0}^{n}A_{k}f(x_{k})$ 中的求积公式用插值公式 $A_{k}=\int_{a}^{b}l_{k}(x)dx$ 表示，则称该数值积分公式为插值型求积公式。

定理：插值型数值积分公式的代数精度至少有n阶精度。
注： $\sum_{k=0}^{n}A_{k}=\sum_{k=0}^{n}\int_{a}^{b}l_{k}(x)dx=\int_{a}^{b}\sum_{k=0}^{n}A_{k}dx=\int_{a}^{b}1dx=b-a$

Newton-Cotes公式

该公式主要研究具体节点的取法，若取[a,b]上n+1个等距节点， $x_{0},x_{1}...x_{n}$ 即第k个节点的值为初始值加上k步步长 $x_{k}=x_{0}+kh$ ，其中 $h=\frac{b-a}{h},x_{0}=a$ ，利用这些节点作n次Lagrange插值多项式，有如下公式:

$\int_{a}^{b}f(x)dx\approx \int_{a}^{b}l_{k}(x)f(x_{k})dx$

$A_{k}=\int_{a}^{b}\frac{(x-x_{0})(x-x_{1})...(x-x_{k-1})(x-x_{k+1})...(x-x_{n})}{(x_{k}-x_{0})(x_{k}-x_{1})...(x_{k}-x_{kk-1})(x_{k}-x_{k+1})...(x_{k}-x_{n})}dx$

令 $x=a+th,x_{k}=a+kh$ ，原式= $\int_{a}^{b}\frac{(b-a)(-1)^n-k}{nk!(n-k)!}\int_{0}^{n}t(t-1)...(t-k+1)(t-k-1)...(t-n)dt$

令 $C_{k}^{n}=\frac{(b-a)(-1)^n-k}{nk!(n-k)!}\int_{0}^{n}t(t-1)...(t-k+1)(t-k-1)...(t-n)dt$ 则为Cotes系数， $A_{k}=(b-a)C_{k}^{n}$

最终Newton-Cotes公式为 $\int_{a}^{b} f(x)dx\approx (b-a)\sum_{k=0}^{n}C_{k}^{n}f(x_{k})$ ，即节点均匀分割的求和公式。

例如当n=1，
$C_{o}^{1}=\frac{(-1)^1}{1!0!1!},\int_{0}^{1}(t-1)dt=\frac{1}{2}$
$C_{1}^{1}=\frac{(-1)^1}{1!1!0!},\int_{0}^{1}tdt=-\frac{1}{2}$
则 $\int_{a}^{b}f(x)dx\approx \frac{b-a}{2}[f(a)+f(b)]$
即在[a,b]区间整体运算，此时计算方法为梯形公式，代数精度为1.

当n=2时，
分别计算 $C_{0}^{2},C_{1}^{2},C_{2}^{2}$ ，
$\int_{a}^{b}f(x)dx\approx \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]$
该公式为Sinpson公式，至少精度为2

n=4时为Cotes公式，具有至少5次代数精度。

定理：n为偶数时，其代数精度至少为n+1次，n表示区间个数，而非节点个数，如5个节点将区间分为四个区间，精度为5。

例题： $f(x)=x^2时，用梯形公式和Simpson计算\int_{0}^{2}x^2dx$ 。
解：
梯形公式 $\int_{0}^{2}x^2dx\approx \frac{2-0}{2}(0+4)=4$
Simpson $\int_{0}^{2}x^2dx\approx \frac{2-0}{6}(0+4+4)=\frac{8}{3}$
可以看到梯形公式的误差还是比较大的，而Simpson甚至得到了精确答案。

性质：
1，对称性。 $C_{k}^{n}=C_{n-k}^n$
2， $\sum_{k=0}^{n}C_{k}^{n}=1$

Newton-Cotes余项：
梯形公式下 $R(f)=\int_{a}^{b}\frac{f''(\xi)}{2}(x-a)(x-b)dx$ ，存在 $\eta \in [0,1]，使R(f)=\frac{1}{2}f''(\eta)\int_{a}^{b}(x-a)(x-b)dx=-\frac{f''(\eta)}{12}(b-a)^3$

同理:
Simpson余项 $R(f)=-\frac{1}{90}(\frac{b-a}{2})^5f^4(\eta)$
Cotes余项 $R(f)=\frac{2(b-a)}{945}(\frac{b-a}{4})^6f^6(\eta)$

复化求积公式

为了克服容格现象，用分段函数代替原函数，用于求积分时同样应用分治思想，划分区间应用梯形或Simpson公式，使用分段线性插值代替原函数，数学公式表示为 $\int_{a}^{b}f(x)dx= \sum_{k=0}^{n-1}\int_{x_k}^{x_{k+1}}f(x)dx\approx \sum_{k=0}^{n-1}I_k$ ，其中 $I_k=\int_{x_k}^{x_{k+1}}f(x)dx\approx \frac{h}{2}[f(x_k)+f(x_{k+1})]$ 或用Simpson

复化的梯形公式

分割区间使用 $I_=\int_{x_k}^{x_{k+1}}f(x)dx\approx \frac{h}{2}[f(x_k)+f(x_{k+1})]$ ， $I_k= \sum_{k=0}^{n-1}\int_{x_k}^{x_{k+1}}f(x)dx \approx 2\frac{h}{2}[f(a)+f(b)]+2\sum_{k=1}^{n-1}f(x_k)$
多个区间分别计算梯形公式后求和

Romberg公式

复化梯形公式方程 $T_n=\frac{h}{2}[f(a)+f(h)+2\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})]$ ，将 $T_n=T_1(h)$ ，由E-M欧拉麦克劳林公式， $T_1(h)-I=\alpha_1h^2+\alpha_2h^4+...+\alpha_kh^2k+...$

后面不写了，反正使用外推法构造出高精度积分公式 $T_{m+1}(h)=\frac{4^m}{4^{m}-1}T_m(\frac{h}{2})-\frac{1}{4^{m}-1}$

计算过程如下图：
Romberg计算过程
例题：利用复化求积公式将[0,1]8等分，计算 $I=\int_{0}^{1}\frac{4+x}{1+x^2}$ 。
解：
复化梯形公式， $I\approx\frac{1}{2×8}[f(0)+f(1)+2\sum_{1}^7f(x_k)]\approx 3,138988$
复化Simpson， $I\approx\frac{1}{6×4}[]f(0)+f(1)+4\sum_{0}^3f(x_{k+\frac{1}{2}})+2\sum_{1}^3f(x_k)]\approx3.141593$

步长越多结果越精确，但应用中具体步长我们难以直接确定，下面将介绍一种能够使步长动态迭代的方法。

变步长的梯形公式

基本思想为：原步长h折半，再利用复化求积公式，不断折半，用前后两次结果做差，如果精度满足要求则停止折半，故该方法又称折半求积公式。

例如： $\int_a^bf(x)dx$ 将[a,b]等分， $h=\frac{b-a}{n},x_k=a+kh$ ，
复化梯形公式 $T_{2n}=\frac{h}{4}\sum_{k=0}^{n-1}[f(x_k)+2f(x_{k+\frac{1}{2}})+f(x_{k+1})]=\frac{T_n}{2}+\frac{h}{2}\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})$ 。
将已复化的梯形公式再折半步长应用复化公式，该公式实际上只需计算 $x_{k+\frac{1}{2}}$ 即新分节点的函数值。

应用中往往使用 $|t_{2k}-T_{2k-1}|<\varepsilon$ ，前后两次变化之差来控制步长。

例题：用变步长的梯形公式计算 $I=\int^1_0\frac{sinx}{x}dx$
解：令 $f(x)=\frac{sinx}{x}$ ，并补充定义 $f (0) = 1$ 使之连续， $T_1=\frac{1}{2}[f(0)+f(1)]\approx0.9207355$
$T_2=\frac{1}{2}T_1+\frac{1}{2}f(\frac{1}{2})\approx0.9397933$

使用该方法的计算量大，收敛慢，为了改进该问题，引入变步长的Romberg公式。

变步长的Romberg公式

其实就是Romberg的折半法，可视化展示如下：
变步长Romberg
公式计算方式如下：

$T_1=\frac{b-a}{2}[f(a)+f(b)]$
$T_2=\frac{T_1}{2}+\frac{b-a}{2}f(\frac{a+b}{2})$
$S_1=\frac{4}{3}T_2-\frac{1}{3}T_1$
$T_4=\frac{T_2}{2}+\frac{b-a}{4}[f(\frac{3a+b}{4})+f(\frac{3b+a}{4})]$
$S_2=\frac{4}{3}T_4-\frac{1}{3}T_2$
$C_1=\frac{16}{15}S_2-\frac{1}{15}S_1$
$T_8=\frac{1}{2}T_4+\frac{b-a}{8}[f(\frac{7a+b}{8})+f(\frac{5a+3b}{8})+f(\frac{5b+3a}{8})+f(\frac{3b+a}{2})]$
$S_4=\frac{4}{3}T_8-\frac{1}{3}T4$
$C_2=\frac{16}{15}S_4-\frac{1}{15}S_2$
$R_1=\frac{64}{63}C_2-\frac{1}{63}C_1$

使用折半+外推法计算Romberg，T使用折半法迭代， $S=\frac{4}{3}T_{k+1}-\frac{1}{3}T_k$ ， $C=\frac{16}{15}S_{k+1}-\frac{1}{15}S_k$ ， $R=\frac{64}{63}C_{k+1}-\frac{1}{63}C_{k}$ 。

例题：使用Romberg公式计算 $\int_2^8\frac{1}{2x}dx$ ，每步计算保留4位小数。
解：
$T_1=\frac{8-2}{2}[f(8)+f(2)]\approx0.9375$
$T_2=\frac{T_1}{2}+\frac{8-2}{2}f(5)\approx0.76875$
$S_1=\frac{4}{3}T_2-\frac{1}{3}T_1=0.7125$
$T_4=\frac{1}{2}T_2+\frac{8-2}{4}[f(3.5)+f(6.5)]=0.7140$
$S_2=\frac{4}{3}T_4+\frac{1}{3}T_2=0.6958$
$C_1=\frac{16}{15}S_2-\frac{1}{15}S_1=0.6947$
$T_8=\frac{1}{2}T_4+\frac{8-2}{8}[f(2.75)+f(4.25)+f(5.75)+f(7.25)]=0.6986$
$S_4=\frac{4}{3}T_8-\frac{1}{3}T_4=0.6934$
$C_2=\frac{16}{15}S_4-\frac{1}{15}S_2=0.6932$
$R_1=\frac{64}{63}C_2-\frac{1}{63}C_1=0.6932$

其他参考资料在此

总结

本章学习了如何利用迭代公式计算积分的近似值：

首先是使用数值积分矩形或梯形来代替原有积分曲线的方法，并用自变量 $x$ 的次数作为代数精度衡量近似的程度。

然后是利用插值多项式求和代替原积分曲线积分，其中点的取法用Newton-Cotes均匀选取。

最后为了计算迭代终点，采用变步长求积动态折半计算，不断逼近真实值，有基本的复化梯形公式和加快迭代收敛速度的Romberg公式，Romberg计算方法最复杂，需要迭代四次计算，但收敛速度快，计算量也没那么大，可以说是计算数值积分最先进的方法。