引言

数值积分的概念

基本思想

数值积分的基本思想是：用多项式函数\(p\)去逼近函数\(f\),然后用近似多项式的积分值去代替原来函数的积分值. \[ \int_{a}^{b}f\,dx\approx\int_{a}^{b}p\,dx \]

<______________________假装是图片_______________________________>

最常用的逼近多项式就是前面所学的插值多项式.

\[ \int_{a}^{b}f\,dx\approx\sum_{i=0}^{n}A_if(x_i)=I_n(f) \] - \(A_i\)为求积系数,它只与节点选取有关,与函数\(f\)无关. - \(I_n(f)\)为数值求积公式.

代数精度

定义: 如求积公式对被积函数\(f(x)=1,x,\cdots,x^m\)都能精确成立,即 \[ \int_{a}^{b}f\,dx=\sum_{i=0}^{n}A_if(x_i)=I_n(f) \] 但对\(f(x)=x^{m+1}\)不能精确成立,即 \[ \int_{a}^{b}x^{m+1}\,dx\neq\sum_{i=0}^{n}A_ix^{m+1}_i \] 则称求积公式有\(m\)次代数精度.

若 \[ I_n(f)=\sum_{i=0}^{n}A_if(x_i) \] 具有\(n\)次代数精度,问:\(A_i=?\) 实际上,分别令\(f=1,x,\cdots,x^n\),有 \[ \int_{a}^{b}f\,dx=\sum_{i=0}^{n}A_if(x_i)=I_n(f) \] 将这\(n+1\)个方程写成线性方程组的形式, \[ \left(\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{0} & x_{1} & \cdots & x_{n} \\ \vdots & \vdots & & \vdots \\ x_{0}^{n} & x_{1}^{n} & \cdots & x_{n}^{n} \end{array}\right)\left(\begin{array}{c} A_{0} \\ A_{1} \\ \vdots \\ A_{n} \end{array}\right)=\left(\begin{array}{c} b-a \\ \frac{b^{2}-a^{2}}{2} \\ \vdots \\ \frac{b^{n+1}-a^{n+1}}{n+1} \end{array}\right) \] 由于这\(n+1\)个节点互异,所以求积系数\(A_i\)是存在且唯一的.通过解上面的线性方程组得到求积系数的方法被称为待定系数法. 结合上面的描述,我们可以很容易得到下面的定理. 定理: 给定\(n+1\)个互异节点,那么总存在相应的求积系数\(A_i\)使得求积公式 \[ \sum_{i=0}^{n}A_if(x_i)=I_n(f) \] 至少具有\(n\)次的代数精度.

插值型求积公式

基本想法: \[ f(x)\approx L_n(x)=\sum_{i=0}^{n}l_{i,n}f(x_i)\\ \int_{a}^{b}f(x)\,dx\approx\int_{a}^{b}\sum_{i=0}^{n}l_{i,n}f(x_i)\,dx=\sum_{i=0}^{n}\int_{a}^{b}l_{i,n}(x)\,dxf(x_i) \]

定义: 若 \[ \sum_{i=0}^{n}A_if(x_i)=I_n(f) \] 其中, \[ A_i=\int_{a}^{b}l_{i,n}(x)\,dx \] 则称\(I_n(x)\)为插值型求积公式.

误差: \[ \begin{split} R_n(f)&=\int_{a}^{b}f(x)\,dx-\int_{a}^{b}l_{i,n}(x)\,dx\\ &=\int_{a}^{b}(f(x)-l_{i,n}(x))\,dx\\ &=\int_{a}^{b}\frac{f^{(n+1)}{(\xi_n)}}{(n+1)!}\pi_{n+1}(x)\,dx \end{split} \]

定理: 求积公式 \[ \sum_{i=0}^{n}A_if(x_i)=I_n(f) \] 为插值型求积公式的充要条件为:\(I_n(x)\)至少具有n次代数精度.

Proof: (必要性) (充分性)

Newton-Cotes 求积公式

节点为等距节点 \[ x_i=x_0+ih,\quad h=\frac{b-a}{n} \] 变形, \[ \begin{split} I_n(f)&=\sum_{i=0}^{n}A_if(x_i)\\ &=(b-a)\sum_{i=0}^{n}\frac{A_i}{b-a}f(x_i)\\ &=(b-a)\sum_{i=0}^{n}C_{i,n}f(x_i)\\ \end{split} \] 由\(x=x_0+th,(0\leq t\leq n)\),可知, \[ \begin{split} C_{i,n}&=\frac{1}{b-a}A_i\\ &=\frac{1}{b-a}\int_{a}^{b}l_{i,n}(x)\,dx\\ &=\frac{1}{b-a}\int_{a}^{b}\prod^{n}_{j=0,i\neq i}\frac{x-x_j}{x_i-x_j}\,dx\\ &=\frac{1}{b-a}\int_{a}^{b}\prod^{n}_{j=0,i\neq i}\frac{t-j}{i-j}h\,dx\\ &=\frac{h}{b-a}\int_{a}^{b}\prod^{n}_{j=0,i\neq i}\frac{t-j}{i-j}\,dx\\ &=\frac{1}{n}\int_{a}^{b}\prod^{n}_{j=0,i\neq i}\frac{t-j}{i-j}\,dx\\ \end{split} \] 而\(C_{i,n}\)被称为Cotes系数(只与\(i,n\)有关) 当\(n=1\)时, \[ C_{0,1}=\frac{1}{2},C_{1,1}=\frac{1}{2}\\ I_1(f)=(b-a)(\frac{f(b)+f(a)}{2}) \] 称为梯形公式. 当\(n=2\)时, \[ C_{0,2}=\frac{1}{6},C_{1,2}=\frac{4}{6},C_{2,2}=\frac{1}{6}\\ I_1(f)=(b-a)(\frac{f(b)+f(a)+4f(\frac{a+b}{2})}{6}) \] 称为Simpson(辛普森)求积公式.

Cotes系数的性质:

• Cotes系数的代数和为1 \[\sum_{i=0}^{n}C_{i,n}=1 \] 由于NC求积公式至少具有\(n\)阶代数精度, \[Let\ f(x)\equiv 1,so\,we\,have\\ \int_{a}^{b}f(x)\,dx=b-a=(b-a)\sum_{i=0}^{n}C_{i,n}f(x_i)=(b-a)\sum_{i=0}^{n}C_{i,n}\]
• 中心对称性质. \[C_{i,n}=C_{n-i,n} \]

Newton-Cotes 公式的误差

Simpson求积公式具有3次代数精度.

定理: Newton-Cotes公式当n为奇数时,有n次代数精度;当n为偶数时,有n+1次代数精度.

Proof:

n为奇数,显然.

n为偶数.\(f(x)=x^{n+1}\) \[ \begin{split} &\int_{a}^{b}f(x)\,dx-(b-a)\sum_{i=0}^{n}C_{i,n}f(x_i)\\ &=\int_{a}^{b}[f(x)-L_n(x)]\,dx\\ &=\int_{a}^{b}\frac{f^{(n+1)}(\xi)}{(n+1)!}\pi_{n+1}(x)\,dx\\ &=\int_{a}^{b}\frac{(n+1)!}{(n+1)!}\pi_{n+1}(x)\,dx\\ &=\int_{a}^{b}\pi_{n+1}(x)\,dx\\ &=\int_{a}^{b}(x-x_0)\cdots(x-x_n)\,dx \end{split} \] 由于是等距节点,所以有,！！！！！！！！！！！！！！！！！！！！ \[ \begin{split} &\int_{a}^{b}\pi_{n+1}(x)\,dx\\ &=\int_{a}^{b}(x-x_0)\cdots(x-x_n)\,dx\\ &=\int_{a}^{b}(x-x_0)\cdots(x-x_n)\,dx \end{split} \]

Newton-Cotes公式的稳定性

定义:(数值稳定性)

定理: 数值求积公式\(I_n(f)=\sum_{i=0}^{n}A_if(x_i)\)为数值稳定的充分必要条件是:\(A_i\geq 0,i=0,\cdots,n\).

Proof

\[ \forall \epsilon>0,\qquad|\sum_{i=0}^{n}A_if(x_i)-\sum_{i=0}^{n}{A_i}\tilde{f}(x_i)|\leq \sum_{i=0}^{n}|A_if(x_i)-\sum_{i=0}^{n}{A_i}\tilde{f}(x_i)|\\ \leq M\sum_{i=0}^{n}{A_i}<\epsilon,as\ long\ as\ M<\delta=\frac{\epsilon}{\sum_{i=0}^{n}{A_i}} \]

复合求积公式

为了避免高阶求积公式的弊端，就要用到插值中的分段思想，即把积分区间分段，再在各个区间上应用低阶的求积公式，称为复合求积公式。

复合梯形求积公式

\[ \begin{split} \int_{a}^{b}f(x)\,dx&=\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx\\ &\approx\sum_{i=0}^{n-1}\frac{x_{i+1}-x_i}{2}(f(x_i)+f(x_{i+1})) \end{split} \]

特别地，考虑等距节点的情况，间距为\(h\), \[ \begin{split} \int_{a}^{b}f(x)\,dx&=\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx\\ &\approx \frac{h}{2}(f(x_0)+2\sum_{i=1}^{n-1}f(x_i)+f(x_n)) \end{split} \]

复合梯形求积公式的误差

上面得出的复合梯形求积公式记为\(T_n(x)\).

\[ \begin{split} \Big|\int_{a}^{b}f(x)\,dx-T_n(x)\Big|&=\Big|\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx-T_n(x)\Big|\\ &=\Big|\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx-\sum_{i=0}^{n-1}\frac{x_{i+1}-x_i}{2}(f(x_i)+f(x_{i+1}))\Big|\\ &\leq \sum_{i=0}^{n-1}\Big|\int_{x_i}^{x_i+1}f(x)\,dx-\frac{x_{i+1}-x_i}{2}(f(x_i)+f(x_{i+1}))\Big|\\ &=\sum_{i=0}^{n-1}\Big|+\frac{(x_{i+1}-x_i)^3}{12}f''(\xi_i)\Big|\\ &=\frac{M}{12}\sum_{i=0}^{n-1}\Big|(x_{i+1}-x_i)^3\Big| \qquad M=\max{f''(\xi_i)}\\ &=\frac{M}{12}h^3\sum_{i=0}^{n-1}(1) \qquad h=\max{(x_{i+1}-x_i)}\\ &=\frac{M}{12}h^2(hn)=\frac{M}{12}h^2(b-a)\qquad if\ nodes\ spaced\ equally \end{split} \]

复合Simpson求积公式

\[ \begin{split} \int_{a}^{b}f(x)\,dx&=\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx\\ &\approx\sum_{i=0}^{n-1}\frac{x_{i+1}-x_i}{6}(f(x_i)+4f(\frac{x_i+x_{i+1}}{2})+f(x_{i+1})) \end{split} \]

特别地，考虑等距节点的情况，间距为\(h\), \[ \begin{split} \int_{a}^{b}f(x)\,dx&=\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx\\ &\approx \frac{h}{6}(f(x_0)+2\sum_{i=1}^{n-1}f(x_i)+f(x_n)+4\sum_{i=1}^{n-1}f(\frac{x_i+x_{i+1}}{2})) \end{split} \]

复合Simpson求积公式的误差

上面得出的复合Simpson求积公式记为\(S_n(x)\).

\[ \begin{split} &\Big|\int_{a}^{b}f(x)\,dx-S_n(x)\Big|\\ &=\Big|\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx-S_n(x)\Big|\\ &=\Big|\sum_{i=0}^{n-1}\int_{x_i}^{x_i+1}f(x)\,dx-\sum_{i=0}^{n-1}\frac{x_{i+1}-x_i}{6}(f(x_i)+4f(\frac{x_i+x_{i+1}}{2})+f(x_{i+1}))\Big|\\ &\leq \sum_{i=0}^{n-1}\Big|\int_{x_i}^{x_i+1}f(x)\,dx-\frac{x_{i+1}-x_i}{6}(f(x_i)+4f(\frac{x_i+x_{i+1}}{2})+f(x_{i+1}))\Big|\\ &=\sum_{i=0}^{n-1}\Big|-\frac{(x_{i+1}-x_i)^5}{2880}f^{(4)}(\xi_i)\Big|\\ &=\frac{M}{2880}\sum_{i=0}^{n-1}\Big|(x_{i+1}-x_i)^5\Big| \qquad M=\max{f^{(4)}(\xi_i)}\\ &=\frac{M}{2880}h^5\sum_{i=0}^{n-1}(1) \qquad h=\max{(x_{i+1}-x_i)}\\ &=\frac{M}{2880}h^4(hn)=\frac{M}{2880}h^4(b-a)\qquad if\ nodes\ spaced\ equally \end{split} \]

Richardson 外插值法和Romberg 算法

上一节得到了复合求积公式，通过外插值方法可以进一步提高复合求积公式的效率，即由低阶求积公式构造高阶求积公式。

Richardson 外插值法

通项公式

\[ \begin{cases} R(j,0)=\Phi(\frac{h}{2^j}),\quad j\geq 0\\ \\ R(j,k)=\frac{2^kR(j,k-1)-R(j,k-1)}{2^{k}-1},\quad j\geq k> 0 \end{cases} \]

Romberg 算法

Euler-Maclaurin展开式

设\(f\in C^{2m}([a,b])\)，那么，

\[ \begin{split} \int_{a}^{b}f(x)\,dx-T_n(f)&=\sum_{r=1}^{n-1}C_rh^{2r}[f^{(2r-1)}(b)-f^{(2r-1)}(a)]+O(h^{2r})\\ &C_r=\frac{}{} \end{split} \]

注意，当\(f\)为周期函数时，

\[ \begin{split} f(b)&=f(a)\\ f^{(2r-1)}(b)&=f^{(2r-1)}(a)\\ \int_{a}^{b}f(x)\,dx-T_n(f)&=O(h^{2r})\\ \end{split} \]

精度已经很高了.

梯形加速公式

Simpson加速公式

Romberg 算法的误差

定理：设\(f\in C^{2m+2}([a,b])\)，那么Romberg 算法的误差估计为：

\[ \int_{a}^{b}f(x)\,dx-T_k^m(f)=O(h^{2m+2}),\qquad k\geq 0,m\geq 0 \]

Romberg 算法的数值稳定性

定理：若Romberg 算法可以展开为

\[ T_k^m(f)=\sum_{i=0}^{2^k}c_{i,m}f(x_i),\qquad k\geq m,m=1,2\cdots \]

则所有系数\(c_{i,m}\)为正数。

Gauss 求积公式

为了简化处理过程,前面都是使用的等距节点作为求积节点,实际上,节点的选取也会影响代数精度.本节将介绍具有更高代数精度的求积方法.

考虑定积分 \[ \int_{a}^{b}\rho(x)f(x)\,dx\approx\sum_{i=0}^{n}A_{i}f(x_i) \]

定理: 求积公式\(\sum_{i=0}^{n}A_{i}f(x_i)\)的代数精度最高不超过\(2n+1\).

Proof: 只要找到一个\(2n+2\)次的多项式求积公式,求积公式不能精确成立. 设 \[ p(x)=\pi_(n+1)^2(x)=\prod_{i=0}^{n}(x-x_i)\in P_{2n+2} \]

\[ \int_{a}^{b}\rho(x)p(x)\,dx=\int_{a}^{b}\rho(x)\prod_{i=0}^{n}(x-x_i)\,dx>0 \]

令一方面,

\[ \sum_{i=0}^{n}A_{i}p(x_i)=0 \]

故,求积公式对\(p(x)\)不能精确成立.

定义: 若由一组节点\(x_i,i=0,\cdots,n\)使得\(\sum_{i=0}^{n}A_{i}f(x_i)\)具有\(2n+1\)次代数精度,则称这组节点为Gauss点,相应的求积公式为Guass求积公式.

常用Gauss求积公式

• Gauss-Legendre 求积公式
• Gauss-Chebyshev 求积公式

Gauss 求积公式的误差估计