数值积分

引入：数值积分的两种思路

思路1：使用Newton-Leibniz公式

\int_a^bf(x){\rm d}x=F(b)-F(a)

问题：一些函数没有初等的原函数

思路2：使用积分的定义

\int_a^bf(x){\rm d}x=\lim\limits_{n\to\infty, h\to0}\sum\limits_{i=0}^n(x_{i+1}-x_i)f(\xi_i)

问题：计算量太大

解决方案：使用机械求积公式

I_n(f)=\sum\limits_{k=0}^n A_kf(x_k)

即：使用 $f(x)$ 的若干个函数值加权求和得到

注为什么这里的记号使用 $I_n(f)$ 而不是 $I_n(x)$ ？因为一般而言求积公式是一种方法，它对任何满足一定条件的函数都是通用的，输入一个函数，输出一个在给定区间上的积分值。积分值与 $f$ 有关，与 $x$ 无关

插值型求积公式

假设对于此函数在区间 $[a,b]$ 内有已知点 $x_0,x_1,\cdots,x_n$ ，其余函数值未知，则根据Lagrange插值公式：

p(x)=\sum\limits_{k=0}^nf(x_k)l_k(x)

于是得到积分公式

I_n(f)=\sum\limits_{k=0}^nf(x_k)\int_a^bl_k(x){\rm d}x

即积分公式中的积分系数

A_k=\int_a^bl_k(x){\rm d}x

这对应着每一个插值点处的脉冲基函数的积分

求积公式的通用讨论

积分余项与代数精度

定义（积分余项） 设使用公式 $I_n(f)$ 计算积分 $I(f)$ ，则

R[f]=I(f)-I_n(f)

称为此积分的积分余项

注它反应了计算的截断误差

若 $I_n(f)$ 为插值函数 $p(x)$ 的积分，则

R[f]=\int_a^b[f(x)-p(x)]{\rm d}x

即此时的积分余项为插值余项的积分，根据插值余项的估计可知

R[f]=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x){\rm d}x

定义（代数精度） 若 $I_n(f)$ 对于被积函数 $f(x)$ 是次数不超过 $m$ 的多项式的情况下都是准确的，对 $m+1$ 次多项式可能不准确，则称该求积公式具有 $m$ 次代数精度

对于梯形公式：

I_1(f)=\frac{b-a}{2}[f(a)+f(b)]

它具有1次代数精度

中矩形公式：

I_0(f)=(b-a)f(\frac{a+b}{2})

也具有1次代数精度，因此代数精度与求积用到的插值点数目并没有必然的关系

定理机械求积公式 $I_n(f)=\sum_{k=0}^nA_kf(x_k)$ 至少有 $m$ 次代数精度等价于此公式对 $f(x)=1,x,\cdots,x^m$ 分别准确

推论求积公式至少 $0$ 次代数精度时，满足 $\sum_{k=0}^nA_k=b-a$

注只要满足上面的等式，求积公式就一定具有0次代数精度，不满足此等式的求积公式可以认为是不可用的，因此上面的等式可以被认为是对可用的求积公式都成立的结论

定理插值型求积公式 $I_n(f)$ 至少具有 $n$ 次代数精度，且至少具有 $n$ 次代数精度的求积公式一定是插值型求积公式

收敛性与稳定性

定义（求积公式的收敛性） 对于机械求积公式 $I_n(f)=\sum_{k=0}^nA_kf(x_k)$ ，设 $h$ 为所有插值点间距的最大值，则求积公式收敛定义为

\lim_{n\to\infty,h\to0}I_n(f)=\int_a^bf(x){\rm d}x

数值积分问题的敏感性分析

设函数从 $f(x)$ 变为 $\tilde f(x)$ ，扰动的大小定义为

\delta=||f(x)-\tilde f(x)||_\infty=\max_{a\leq x\leq b}|f(x)-\tilde f(x)|

则结果误差

\left|\int_a^bf(x){\rm d}x-\int_a^b\tilde f(x){\rm d}x\right|\leq\int_a^b\left|f(x)-\tilde f(x)\right|{\rm d}x\leq(b-a)\delta

于是问题的绝对条件数的上限为 $(b-a)$

机械求积公式的敏感性分析

当 $f(x)$ 变为 $\tilde f(x)$ 时，

|I_n(f)-I_n(\tilde f)|=\left|\sum\limits_{k=0}^nA_k[f(x_k)-\tilde f(x_k)]\right|

\leq\sum\limits_{k=0}^n|A_k|\cdot|f(x_k)-\tilde f(x_k)|\leq \left(\sum\limits_{k=0}^n|A_k|\right)\epsilon

其中

\epsilon = \max_{0\leq k\leq n}|f(x_k)-\tilde f(x_k)|\leq \delta

由于我们已知

\sum\limits_{k=0}^nA_k=b-a

因此若所有 $A_k$ 均正，我们有

|I_n(f)-I_n(\tilde f)|\leq(b-a)\delta

否则无法得出有效的放缩，因此

定义（机械求积公式的稳定性） 若所有 $A_k$ 均正，则称机械求积公式 $I_n(f)$ 是稳定的

Newton-Cotes公式

使用插值型求积公式，让节点 $x_k$ 在区间上均匀分布，即 $h=(b-a)/n$ ， $x_k=a+kh$

此时得到的求积公式为

I_n(f)=\sum\limits_{k=0}^nA_kf(x_k)

A_k=\int_a^b\frac{(x-x_0)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)}{(x_k-x_0)\cdots(x_k-x_{k-1})(x_k-x_{k+1})\cdots(x_k-x_n)}{\rm d}x

这称为 $n$ 阶Newton-Cotes公式

此时作换元 $x=a+th$ ，则

A_k=\int_0^n\prod_{j=0,j\neq k}^n\left(\frac{t-j}{k-j}\right){\rm d}(a+th)=\int_0^n\prod_{j=0,j\neq k}^n\left(\frac{t-j}{k-j}\right)\frac{b-a}{h}{\rm d}t

=(b-a)C_k^{(n)}

其中 $C_k^{(n)}$ 与积分区间无关，称为 $n$ 阶N-C公式的Cotes系数，其满足

\sum\limits_{k=0}^nC_k^{(n)}=1

C_k^{(n)}=C_{n-k}^{(n)}

当 $n=8$ 时，存在负的 $C_k^{(n)}$ ，此时N-C公式不稳定

三种常用的低阶N-C公式：

n=1:\quad T(f)=(b-a)\left[\frac12f(a)+\frac12f(b)\right]

n=2:\quad S(f)=(b-a)\left[\frac16f(a)+\frac46f(\frac{a+b}2)+\frac16f(b)\right]

n=4:\quad C(f)=(b-a)\left[\frac7{90}f(x_0)+\frac{32}{90}f(x_1)+\frac{12}{90}f(x_2)+\frac{32}{90}f(x_3)+\frac{7}{90}f(x_4)\right]

$n=1,2,4$ 时分别称为梯形公式、Simpson公式、Cotes公式

定理（N-C公式的额外代数精度） $n$ 为偶数时， $n$ 阶N-C公式至少有 $n+1$ 阶代数精度

注这说明偶数阶N-C公式相比一般的插值型积分公式有一阶额外的代数精度

证明只需要讨论 $f(x)=x^{n+1}$ 时的积分余项

R[f]=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x){\rm d}x=\int_a^b\omega_{n+1}(x){\rm d}x=h^{n+2}\int_0^n\prod\limits_{j=0}^n(t-j){\rm d}t

此时考虑积分

\int_0^n\prod\limits_{j=0}^n(t-j){\rm d}t

作换元 $t=u+n/2$ ，则

\int_{-n/2}^{n/2}H(u){\rm d}u

其中

H(u)=\prod_{j=0}^n(u+n/2-j)=\prod_{j=-n/2}^{n/2}(u-j)

此时 $H(u)$ 是奇函数，因此上述积分为 $0$ ，故偶数阶N-C公式的对 $x^{n+1}$ 也准确

注这说明偶数阶N-C公式的额外代数精度来自于其插值点选取的对称性

注因此最常用的N-C公式为 $n=1,2,4$ 的情况， $n=3$ 的情形一般不用，它与 $n=2$ 的代数精度相同但计算量更大

低阶N-C公式的积分余项

梯形公式

R_T=I(f)-T(f)=\int_a^b\frac{f''(\xi)}{2}(x-a)(x-b){\rm d}x

此时应用第一积分中值定理可知

R_T=\frac{f''(\eta)}{2}\int_a^b(x-a)(x-b){\rm d}x=-\frac{f''(\eta)}{12}(b-a)^3

Simpson公式

R_S=-\frac{f^{(4)}(\eta)}{2880}(b-a)^5

N-C公式的稳定性和收敛性

稳定性

$n=8$ 的N-C公式不稳定，且 $n\geq10$ 时必定不稳定

收敛性

由于龙格现象，即高次多项式插值不能很好地收敛于原始函数，N-C公式在 $n$ 增大的情况下并不一定收敛到准确解

结论

一般只使用 $n=1,2,4,6$ 的N-C公式

主要计算量为 $n+1$ 次函数求值，其余计算量都可以通过查Cotes系数表省去

复合求积公式

高阶N-C公式不稳定，于是并非阶数越高越好

在函数插值中，同样也并非越高次插值越好，可以采用分段低次插值

于是，可以从分段插值导出复合求积公式

复合梯形公式

将 $[a,b]$ 等分为 $n$ 个子区间，对每个子区间采用梯形公式，则

T_n(f)=\sum\limits_{k=0}^{n-1}\frac h2[f(x_k)+f(x_{k+1})]

其中 $h=(b-a)/n$

其积分误差可以估计为：

\sum\limits_{k=0}^{n-1}\frac{f''(\eta_k)}{12}\left(\frac{b-a}{n}\right)^3\approx\frac{f''(\eta)(b-a)^3}{12n^2}

此截断误差随区间长度 $h$ 减小按 $O(h^2)$ 减小，这被称为具有2阶准确度

类似可以定义 $p$ 阶准确度，定义略

上面的求积公式也是一种机械求积公式：

T_n(f)=\frac h2\left[f(a)+2\sum\limits_{k=0}^{n-1}f(x_k)+f(b)\right]

没有负数系数，故此求积公式稳定且收敛

复合Simpson公式

S_n=\frac h6\sum\limits_{k=0}^{n-1}[f(x_k)+4f(x_{k+1/2})+f(x_{k+1})]

=\frac h6\left[f(a)+f(b)+4\sum\limits_{k=0}^{n-1}f(x_{k+1/2})+2\sum\limits_{k=1}^{n-1}f(x_k)\right]

其中 $f(x_{k+1/2})$ 为 $f(x_k)$ 和 $f(x_{k+1})$ 的中点

$S_n$ 具有稳定性、收敛性、3次代数精度、4阶准确度

注复合求积公式的代数精度一般和对应的简单求积公式相同

注代数精度描述的是对怎样的函数能够无误差，精确度描述的是对任意的函数能够有多小的误差；代数精度是静态的代数结论，精确度是随区间逐渐细分的渐进动态过程

步长折半策略

在实际计算积分时，我们往往有一个给定的精度需求，即可容忍的最大误差，但误差与步长 $h$ 只有渐进表达式，在实际计算中往往动态确定步长

此时可以采用步长折半策略：每次计算后，如果没有达到给定的精度要求，则将 $h$ 变为上一次的一半继续计算，此时有一半的插值点已经完成计算，只需要重新计算

具体地，折半前：

T_n=\frac h2\left[f(a)+2\sum\limits_{k=1}^{n-1}f(x_k)+f(b)\right]

折半后：

T_{2n}=\frac h4\left[f(a)+2\sum\limits_{k=1}^{2n-1}f(x_{k/2})+f(b)\right]

=\frac12T_n+\frac h2\sum\limits_{k=0}^{n-1}f(x_{k+1/2})

即：只需计算折半时新增节点的函数值，且这种重用不需要存储之前计算好的函数值

注复合Simpson公式计算时重用效率会低一些，但仍然可以重用；而中矩形公式一般不用于复合求积，就是因为复合中矩形公式无法在步长折半时复用计算

理查森外推方法

回忆之前得到的复合梯形公式的余项：

I(f)-T_n=-\frac{h^2(b-a)}{12}f''(\eta)=O(h^2)

事实上，这个 $O(h^2)$ 是某个对每个 $f$ 确定的 $h$ 的函数

定理设被积函数 $f(x)\in C^\infty[a,b]$ ，记 $T(h)$ 为积分步长为 $h$ 时复合梯形公式的结果，则

T(h)=I(f)+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots

其中系数 $\alpha_l$ 与 $h$ 无关

注于是，我们可以通过一些 $h$ 比较大的近似计算的结果来反推出 $\alpha_l$ 的近似值，从而利用上面的余项展开式来得到更精确的结果，这就是理查森外推法的思想

然而， $\alpha_l$ 虽然与 $h$ 无关，但与 $f$ 有关，因此其精确值一般情况下也是不可知的

我们可以在步长折半的过程中通过外推法得到更高阶的求积公式，例如，设折半前步长为 $h$ ，计算结果为

T_0^{(n)}=I+\alpha_1h^2+\alpha_2h^4+\cdots

折半后的结果为

T_0^{(n+1)}=I+\frac{\alpha_1}4h^2+\frac{\alpha_2}{16}h^4+\cdots

使用这两个公式拼凑，可以抵消掉 $h^2$ 项，得到

T_1^{(n)}=\frac{4T_0^{(n+1)}-T_0^{(n)}}{3}=I+O(h^4)

于是得到了一个更高阶（此处指的准确度而非代数精度）的求积公式

在实际计算中，可以列三角形数表 $T_k^{(n)}$ 计算

注三角形数表上下标的含义

$T_k^{(n)}$ 由 $T_{k-1}^{(n)}$ 和 $T_{k-1}^{(n+1)}$ 外推得到

$k$ 表示得到这个公式的外推阶数，对应 $2k+2$ 阶准确度； $n$ 表示得到这个公式的所有基础公式的最小二分次数，在这个公式中的 $h$ 相等于 $(a-b)/2^n$

理查森外推算法的问题：

对函数有光滑性要求
外推一次牺牲了步长，例如 $T_1^{(0)}$ 需要使用步长为 $(a-b)/2$ 的点来计算，但其对应的余项公式中的 $h$ 实际为 $a-b$

使用理查森外推的积分公式也叫Romberg积分公式

高斯求积公式

Newton-Cotes公式中最重要的前提是选取等距插值节点，但这其实是不必要的。N-C公式中，我们通过对称选取插值节点，对偶数阶公式取得了更高的代数精度，但事实上，选取恰当的插值节点还能得到更高的代数精度

考虑机械求积公式

I_n=\sum\limits_{k=0}^nA_kf(x_k)

此公式中积分系数和积分节点一共有 $2n+2$ 个待定参数，理论上可以达到至少 $2n+1$ 次的代数精度

定义满足上面说的最高代数精度的求积公式称为（ $n$ 阶）高斯求积公式，此时的积分节点称为高斯点

由于需要满足 $2n+1$ 阶代数精度的方程是通过指定被积函数为 $1,x,\cdots,x^{2n+1}$ 得到的，因此这些待定参数与被积函数 $f$ 无关，可以实现求出；但由于待定参数中包含高斯点，高斯点与积分上下限有关

例（1阶高斯公式）

假设积分限为 $[-1,1]$ ，所求积分公式为

G_1(f)=A_0f(x_0)+A_1f(x_1)

将 $f=1,x,x^2,x^3$ 分别代入：

A_0+A_1=2\\ A_0x_0+A_1x_1=0\\ A_0x_0^2+A_1x_1^2=\frac23\\ A_0x_0^3+A_1x_1^3=0

解得

G_1(f)=f(-\frac{\sqrt 3}3)+f(\frac{\sqrt 3}3)

一般高斯求积公式求法

高斯求积公式是插值型求积公式，因此只要能确定高斯点 $x_k$ 就可以确定系数 $A_k$ ：

A_k=\int_a^bl_k(x){\rm d}x

主要考虑如何求解高斯点

考虑带权积分：

I=\int_a^bf(x)\rho(x){\rm d}x

这里的权函数 $\rho(x)$ 满足在积分限内非负、可积

上述所有积分理论都可以应用到带权积分中，对每个权函数可以求出一套高斯点，下面讨论对带权积分的通用求法

定理（高斯点的充要条件） $x_k,\ k=0,1,\cdots,n$ 是高斯点的充要条件是多项式

\omega_{n+1}(x)=(x-x_0)(x-x_1)\cdots(x-x_n)

与 $\forall P(x)\in\mathbb P_n$ 正交，即

\int_a^bP(x)\omega_{n+1}(x)\rho(x){\rm d}x=0

注上面的正交指的是，在给定积分限和权函数的积分体系下正交

证明 - 必要性

当 $x_k$ 是高斯点时，高斯积分有 $2n+1$ 次代数精度，故对任何不超过 $n$ 次的多项式 $P$ ， $P(x)\omega_{n+1}(x)$ 不超过 $2n+1$ 次，故

I_n(\omega_{n+1}P)=I(\omega_{n+1}P)

左侧等于

\sum\limits_{k=0}^nA_kP(x_k)\omega_{n+1}(x_k)=0

故右侧也为 $0$

证明 - 充分性

对于任何 $f(x)\in\mathbb P_{2n+1}$ ，可以将其除以多项式 $\omega_{n+1}$ ，得到

f(x)=P_n(x)\omega_{n+1}(x)+Q_n(x)

其中

P_n,Q_n\in\mathbb P_n

于是

\int_a^bf(x)\rho(x){\rm d}x=\int_a^bP_n(x)\omega_{n+1}(x)\rho(x){\rm d}x+\int_a^bQ_n(x)\rho(x){\rm d}x

右侧第一项根据正交性条件可知为零；而由于 $I_n$ 是 $x_0,\cdots,x_n$ 确定的插值型求积公式，至少有 $n$ 阶代数精度，故有

\int_a^bf(x)\rho(x){\rm d}x=\int_a^bQ_n(x)\rho(x){\rm d}x=\sum\limits_{k=0}^nA_kQ_n(x_k)

由于 $\omega_{n+1}$ 在所有插值点处取 $0$ ，可知 $f(x_k)=Q_n(x_k),\forall k$ ，于是

\int_a^bf(x)\rho(x){\rm d}x=\sum\limits_{k=0}^nA_kf(x_k)=I_n

结论高斯点就是 $n+1$ 次正交多项式的零点

注正交多项式与权函数有关；根据正交多项式的性质，正交多项式的零点都是积分区间内的单实根

高斯求积公式具有稳定性，这是因为：

A_k=\sum\limits_{j=0}^nA_jl_k^2(x_j)=\int_a^bl_k^2(x)\rho(x){\rm d}x\geq0

注上面第一个等号是由于 $l_k$ 只在对应插值点取非零值而得到的天然相等；第二个等号是由于高斯公式具有 $2n+1$ 阶代数精度，而 $l_k^2$ 是 $2n$ 阶多项式；注意对一般的插值求积公式，第二个等号未必满足

同时，高斯求积公式还具有收敛性，其余项为

R_n[f]=\frac{f^{(2n+2)}(\eta)}{(2n+2)!}\int_a^b\omega_{n+1}^2(x)\rho(x){\rm d}x

当积分区间为 $[-1,1]$ ，权函数为 $\rho(x)=1$ 时，正交多项式退化为勒让德多项式，此时的高斯点即为勒让德多项式的零点，此时的高斯积分公式称为高斯-勒让德公式

此公式的 $x_k$ 和 $A_k$ 可以通过查表得到，比较有实用价值

数值微分

数值微分使用有限差分近似计算导数，常用的差分公式有：

向前差分

D_f(h)=\frac{f(x+h)-f(x)}{h}

向后差分

D_b(h)=\frac{f(x)-f(x-h)}{h}

中心差分

D_c(h)=\frac{f(x+h)-f(x-h)}{2h}

注上述三个公式分别相当于 $[x,x+h]$ 、 $[x-h,x]$ 、 $[x-h,x+h]$ 上的一阶差商，由于差商与导数有关系
$f[x_0,\cdots,x_n]=\frac{f^{(n)}(\xi)}{n!}$
对一阶差商的情形就是 $f'(\xi)$ ，其中 $\xi$ 属于上述对应的区间，这就是差分近似的差商解释

高阶差分近似也可以用差商解释

由于

f(x+h)=f'(x)h+\frac{f''(x)}2h^2+\frac{f'''(\xi_1)}6h^3

f(x-h)=-f'(x)h+\frac{f''(x)}2h^2-\frac{f'''(\xi_2)}{6}h^3

可知

D_c(h)=f'(x)+\frac{f'''(\xi)}6h^2=f'(x)+O(h^2)

向前差分和向后差分在渐进准确度上都不如中心差分

中心差分的误差分析

\epsilon_{\rm tot}=\frac{Mh^2}6+\frac\epsilon h

其中第一项为截断误差， $M$ 为 $|f'''(\xi)|$ 的上界；第二项为传递误差， $\epsilon$ 为计算 $f(x)$ 的误差限

随着 $h$ 变小，第一项减小，第二项增大，最小值在

h=\sqrt[3]{\frac{3\epsilon}M}

时取到，此时可以求得最准确的解

二阶中心差分

f''(x)\approx G_c(h)=2f[x-h, x, x+h]=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}

准确度为2阶

插值型求导公式

类似插值型求积公式，插值型求导公式的思路为利用插值函数的导数近似原始函数的导数

采用Lagrange插值：

L_n(x)=\sum\limits_{k=0}^nf(x_k)l_k(x)\approx f(x)

则导数近似为

f^{(i)}(x)\approx \sum\limits_{k=0}^nf(x_k)l_k^{(i)}(x)

增加插值节点可以得到更高准确度的导数公式，或是求更高阶导数的公式；在给定插值点个数的情况下，准确度和导数阶数提高一个一般会降低另一个

数值微分的外推算法

类似数值积分，中心差分公式也存在某种误差展开式：

D_c(h)=f'(x)+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_kh^{2k}+\cdots

其中 $\alpha_k$ 与 $h$ 无关

此时可以类似积分中做的方法，使用理查森外推构造更准确的公式

这里本质上也是在通过提高插值点个数来提高精度