工程数学 计算方法 第五章 数值积分

数值积分

数值微分

微分

定义：
$$f\left(x\right)=\underset{h\to 0}{\lim }\frac{f\left(a+h\right)-f\left(h\right)}{h}$$
实际运算中会遇到的困难：

1. f(x)以函数表的形式给出
2. f(x)复杂，不便求导
3. 计算机不会算极限（编程/计算方法）

数值微分：近似求导

只知函数表，计算微分：用函数数值的线性组合来表示微分 → 数值微分。
即找Ai使得
$$\sum _{i=0}^N{A}_{i}f(x_i)\approx f^{\prime }(a)$$

差商法

原理

斜率近似导数。

已知A(a,f(a))，B(a+h,f(a+h))，C(a-h,f(a-h))，求A点处的导数值:用斜率近似导数.
AB的斜率：向前差商
$$f^{\prime }(x)=\frac{f(a+h)-f(a)}{h}=f[a,a+h]$$

AB的斜率：向后差商
$$f^{\prime }(x)=\frac{f(a)-f(a-h)}{h}=f[a-h,a]$$
BC的斜率：中心差商
$$f^{\prime }(x)=\frac{f(a+h)-f(a-h)}{2h}=G(h)$$

误差分析

方法：泰勒公式展开
$$f(a+h)=f(a)+h{f}^{\prime }(a)+\frac{h^2}{2!}{f}^{\prime \prime }(\xi )$$
向前差商：
$$R(x)=f^{\prime }(a)-\frac{f(a+h)-f(a)}{h}=-\frac{h}{2}{f}^{\prime \prime }(\xi )=O(h)$$
向后差商：
$$R(x)=f^{\prime }(a)-\frac{f(a)-f(a-h)}{h}=-\frac{h}{2}{f}^{\prime \prime }(\xi )=O(h)$$
中心差商：
$$\begin{gathered} R(x)=f^{\prime }(a)-G(a)=\frac{h^2}{12}(f^{\prime \prime \prime }({\xi }_1)+f^{\prime \prime \prime }({\xi }_2))=\frac{h^2}{6}{f}^{\prime \prime \prime }(\xi )=O(h^2) \\ \xi ,{\xi }_1,{\xi }_2\in (a-h,a+h) \end{gathered}$$

1. 由截断误差，步长h越小，精度越高
2. 但h越小，f(a+h)与f(a-h)越接近，由于舍入误差，可能会导致结果为0的情况出现

优化：变步长算法

二分步长，误差采用事后估计法自动选择步长
记
$$D_1=G(h),D_2=G(\frac{h}{2})$$
则
$$\begin{gathered} f^{\prime }(a)-D_1\approx O(h^2),f^{\prime }(a)-D_2\approx O(\frac{h}{2}{)}^2 \\ \Rightarrow \frac{f^{\prime }(a)-D_2}{f^{\prime }(a)-D_1}\approx \frac{1}{4} \\ \Rightarrow f^{\prime }(a)-D_2\approx \frac{1}{3}(D_2-D_1) \\ （|D_2-D_1|⩽\epsilon 时达到精度要求） \end{gathered}$$
事后误差估计法：
在实际计算过程中，f’(x)是不可知的，f’(a)也就是不可知的，故不可用f’(a)-G(a)来计算误差。这里采用计算结果来估计误差，如果误差的估算值小于等于要求,即|D1-D2|≤ε时，便认为已经达到了精度要求。
先采用起始h带入计算，然后对h二分，每次取半。当精度达到要求时停止计算。

外推：提升精度
$$\begin{gathered} f^{\prime }(a)-D_2\approx \frac{1}{3}(D_2-D_1) \\ \Rightarrow f^{\prime }(a)\approx \frac{4}{3}{D}_2-\frac{1}{3}{D}_1=G_1 \\ {G}_1(h)=\frac{4}{3}G(\frac{h}{2})-\frac{1}{3}G(h) \\ 精度：f^{\prime }(a)-G_1(h)\approx (h^4) \end{gathered}$$
还可以进一步外推：
$$\begin{gathered} G_2(h)=\frac{16}{15}{G}_1(\frac{h}{2})-\frac{1}{15}{G}_1(h) \\ {G}_n(h)=\frac{4^n}{4^n-1}{G}_{n-1}(\frac{h}{2})-\frac{1}{4^n-1}{G}_{n-1}(h) \end{gathered}$$

插值法

除了用几何意义近似，求微分也可以知函数表达式套公式→插值获得函数表达式。

已知f(x)函数表，构造f(x)的Lagrange插值多项式Ln(x)
$$\begin{gathered} f(x)\approx L_n(x) \\ {R}_n(x)=\frac{f^{(n+1)}({\xi }_x)}{(n+1)!}\prod _{i=0}^n(x-x_i) \end{gathered}$$
得近似求导公式：
$$f^{(k)}(x)\approx L_n^{(k)}(x)$$
当k=1时：
$$\begin{gathered} f^{\prime }(x)\approx L_n^{\prime }(x) \\ 误差项：R(x)=\frac{d}{dx}(\frac{f^{(n+1)}({\xi }_x)}{(n+1)!}\prod _{i=0}^n(x-x_i)) \end{gathered}$$

当x=xj时：（xj为函数节点。积形式求导为各项依次求导乘余项求和，x=xj时含xj项全部为0，只余下对（x-xj）求导的那一项）

$$\begin{gathered} f^{\prime }(x_j)\approx \sum _{i=0}^n{l}_i^{\prime }(x_j)f(x_i),j=0,1,...,n \\ 误差项：R(x_j)=\prod _{\begin{array}{c}i=0\\ i\ne j\end{array}}^n(x_j-x_i)\frac{f^{(n+1)}({\xi }_x)}{(n+1)!} \end{gathered}$$

注意：无论是差商法还是插值法，都是只有函数节点值求函数的微分，其计算过程都是对函数节点值的线性组合。

 

数值积分

机械求积公式

$$\begin{gathered} I={\int }_a^{b}f(x)dx \\ 估算：I\approx I_n=\sum _{i=0}^n{A}_{i}f(x_i) \\ 误差/余项：R_n(x)=I-\sum _{i=0}^n{A}_{i}f(x_i) \end{gathered}$$

1. 没有f(x)的具体表达式求f(x)的积分，就是求已知函数节点值的线性组合。
2. 求积系数Ak的选取仅与节点xi的选取有关，与f(x)无关：机械。

只学过给函数具体形式的积分
→把函数表变成函数表达式
→插值

插值型数值积分

在[a,b]上去n+1点做n次插值多项式：
$$\begin{gathered} L_n(x)=\sum _{k=0}^{n}f(x_k)l_k(x) \\ {\int }_a^{b}f(x)dx\approx \sum _{k=0}^{n}f(x_k){\int }_a^b{l}_k(x)dx\triangleq \sum _{k=0}^{n}f(x_k)A_k \\ {A}_k={\int }_a^b{l}_k(x)dx={\int }_a^b\prod _{j\ne k}\frac{x-x_j}{x_k-x_j}dx \end{gathered}$$

代数精度⭐

对于插值型数值积分：Ak使插值结果精确成立的条件是什么？精确与否如何衡量？误差如何衡量？

衡量误差：
$$|I-D|⩽\epsilon$$
但实际上无法得知精确值，也没有近似，所以选择换一种衡量尺度：阶数。

定义：

若某求积公式的误差$R[f]$满足：$R[x^k]=0$对任意$k⩽n$阶多项式成立，且$R[x^{n+1}]\ne 0$对某个n+1阶多项式成立，则称其代数精度为n。

求代数精度：
依次带入$f=1,x,x^2,x^3,\cdots ,x^n,\cdots$，最后一个使$I\approx I_n={\sum }_{i=0}^n{A}_{i}f(x_i)$精确成立的幂次即为代数精度。精确成立即要求带入$f$后左右两边形式完全相等。

代数精度高意味着$I=I_n$更大的成立范围。代数精度n时对n阶及以下的多项式时精确的。n越大，可精确求得的范围越大。
$$形如\sum _{k=0}^n{A}_{k}f(x_k)的求积公式至少有n次代数精度\iff A_k={\int }_a^b{l}_k(x)dx$$

节点等距分布

当节点等距分布时插值结果代数精度可能大于n。
$$\begin{gathered} x_i=a+ih,h=\frac{b-a}{n},i=0,1,...,n \\ {A}_i={\int }_{x_0}^{x_n}\prod _{j\ne i}\frac{x-x_i}{x_i-x_j}dx,令x=a+th \\ {A}_i={\int }_0^n\prod _{j\ne i}\frac{(t-j)h}{(i-j)h}hdt=(b-a)\frac{(-1{)}^{n-i}}{ni!(n-i)!}{\int }_0^n\prod _{i\ne j}(t-j)dt=(b-a)C_i^{(n)} \\ 其中C_i^{(n)}=\frac{(-1{)}^{n-i}}{ni!(n-i)!}{\int }_0^n\prod _{i\ne j}(t-j)为Cotes系数，可直接查表 \\ 故等距分布时：{\int }_a^{b}f(x)dx\approx (b-a)\sum _{n=0}^N{B}_{n}f(a+nh) \end{gathered}$$
n为奇数时代数精度为n，n为偶数时代数精度为n+1。

与龙格现象类似，这个方法在n≤6时近似相等，在n≥7时误差较大。

复合求积

高次插值有Runge现象：分段低次插值。

复合求积基础方法

复合梯形公式：每次选择两个点来插值计算，即计算梯形的面积。
$$\begin{gathered} {\int }_a^{b}f(x)dx\approx \frac{b-a}{2}[f(a)+f(b)] \\ {\int }_a^{b}f(x)dx\approx \frac{h}{2}\left[f(a)+2\sum _{n=1}^{N-1}f(x_n)+f(b)\right] \end{gathered}$$
复化Simpson公式：每次选择三个点来插值计算，即曲边梯形曲边为抛物线。
$$\begin{gathered} {\int }_a^{b}f(x)dx\approx \frac{b-a}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right] \\ {\int }_a^{b}f(x)dx\approx \frac{h}{6}\left[f(a)+4\sum _{n=0}^{N-1}f(x_{x+\frac{1}{2}})+2\sum _{n=1}^{N-1}f(x_n)+f(b)\right] \end{gathered}$$
为方便编程，记${\tilde{x}}_{2n}=x_n,{\tilde{x}}_{2n-1}=x_{x-\frac{1}{2}},n=1,2,...,N$则
$$\begin{gathered} {\int }_a^{b}f(x)dx\approx \frac{\tilde{h}}{3}\left\{f(a)-f(b)+2\sum _{n=1}^N[2f({\tilde{x}}_{2n-1})+f({\tilde{x}}_{2n})]\right\} \\ 其中\tilde{h}=\frac{b-a}{2N},{\tilde{x}}_n=a+n\tilde{h},n=1,2,...,2N \end{gathered}$$

收敛速度与误差估计

与代数精度类似，选择阶数来衡量。此处衡量的是误差，类似高阶无穷小的形式。
$$p阶收敛\iff \underset{h\to 0}{\lim }\frac{R[f]}{h^p}=C<\infty$$
收敛阶数越高，同计算量下精度更高。

给定精度需求，如何取n(采样点数)？
n*h=区间长度，所以也就是确定步长h的问题。
若h太大，则结果不精确不容易满足精度要求。若h太小，则计算量太大。
如何既满足精度要求又使得计算量可以接受？
类似变步长算法：从二分开始不断对分最终使得结果满足精度需求：Romberg积分。

龙贝格积分

1. 确定基础步长
2. 二分直至满足精度要求

问题：

• 如何得知达到了精度要求？(精确值I未知)
• 计算量能否优化？

误差的事后估计法

类似微分中的事后估计法。
$$\begin{gathered} I-T_n=O(h^2),I-T_{2n}=O(\frac{h}{2}{)}^2 \\ \therefore \frac{I-T_{2n}}{I-T_n}=\frac{O(\frac{h}{2}{)}^2}{O(h^2)}\approx \frac{1}{4} \\ \therefore I-T_{2n}\approx \frac{1}{3}(T_{2n}-T_n) \\ 若|T_{2n}-T_n|<\epsilon ,则|I-T_{2n}|<\epsilon \end{gathered}$$

简化计算：递推

感觉记着就完了：
梯形公式：
$$T_{2n}=\frac{1}{2}{T}_n+\frac{h}{2}\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}}),其中h=\frac{b-a}{2^{k-1}}$$
注意公式中的h为上一步的步长，即折半前的步长。
每次步长折半即进行一次修正。修正值为 已知值取半加上新采样点的h/2倍。

外推：
$$\begin{gathered} I-T_{2n}\approx \frac{1}{3}(T_{2n}-T_n) \\ \therefore I\approx \frac{4}{3}{T}_{2n}-\frac{1}{3}{T}_n=S_n \end{gathered}$$
即复化Simpson公式。
未增加计算量提高了计算精度。
注意虽然外推得的公式是I的近似值，但与Simpson公式是精确相等的，二者都为I的近似公式。

余下外推类似：
$$\begin{gathered} \therefore \frac{I-T_{2n}}{I-T_n}=\frac{O(\frac{h}{2}{)}^4}{O(h^4)}\approx \frac{1}{16} \\ \therefore I-T_{2n}\approx \frac{1}{15}(T_{2n}-T_n) \\ \therefore I\approx \frac{16}{15}{T}_{2n}-\frac{1}{15}{T}_n=C_n \end{gathered}$$
即复化Cotes公式。

用类似方法可以不断外推提高精度。

整合：

⭐Romberg积分的步骤⭐

先增加分点，然后外推提高精度。
计算步骤：

1. $$T_1=\frac{(b-a)}{2}[f(b)+f(a)]$$

2. $$T_2=\frac{T_1}{2}+\frac{b-a}{2}f(\frac{a+b}{2})$$

3. $$S_1=\frac{4T_2-T_1}{4-1}$$

4. $$T_4=\frac{T_2}{2}+\frac{\frac{b-a}{2}}{2}\sum _{k=0}^{1}f(a+\frac{\frac{b-a}{2}}{2}+\frac{b-a}{2}k)$$

5. $$S_2=\frac{4T_4-T_2}{4-1}$$

6. $$C_1=\frac{4^2{S}_2-S_1}{4^2-1}$$

7. $$T_8=\frac{T_4}{2}+\frac{\frac{\frac{b-a}{2}}{2}}{2}\sum _{k=0}^{3}f(a+\frac{\frac{\frac{b-a}{2}}{2}}{2}+\frac{\frac{b-a}{2}}{2}k)$$

8. $$S_4=\frac{4T_8-T_4}{4-1}$$

9. $$C_2=\frac{4^2{S}_4-S_2}{4^2-1}$$

10. $$R_1=\frac{4^3{S}_2-C_1}{4^3-1}$$

11. $$.......$$

判断精度：
$$|Y_1-X_1|<\epsilon$$

1、3、6、10中有满足精度要求的即可停止计算。
每外推一次收敛阶提高2阶。

书上的写法：
$$\begin{gathered} \{\begin{array}{rc}{T}_0^{(k)}=\frac{1}{2}{T}_0^{(k-1)}+\frac{b-a}{2^k}{\sum }_{n=1}^{2^{k-1}}f\left[a+\frac{2n-1}{2^k}(b-a)\right],k=1,2,...& \\ & \\ T_m^{(l)}=\frac{4^m{T}_{m-1}^{(l+1)}-T_{m-1}^{(l)}}{4^m-1},l=0,1,...,m=1,2,...& \end{array} \\ 终止准则：|T_m^{(0)}-T_{m-1}^{(0)}|<\epsilon ，取积分逼近值为T_m^{(0)} \end{gathered}$$

应试

Romberg积分的步骤
顺序就是算法，顺序不对就是错误。

 

代数精度（选择题）