数值微分与数值积分

1 数值微分¶

1.1 差商公式¶

数值微分的目标是：已知函数 \(f(x)\) 在若干节点上的函数值，近似计算导数值。

由导数定义

\[ f'(x)=\lim_{h\to0}\frac{f(x+h)-f(x)}{h} \]

可以取有限但较小的 \(h\)，得到常见差商公式。

向前差商¶

\[ f'(x)\approx \frac{f(x+h)-f(x)}{h} \]

向后差商¶

\[ f'(x)\approx \frac{f(x)-f(x-h)}{h} \]

中心差商¶

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

1.2 截断误差¶

由 Taylor 展开

\[ f(x+h)=f(x)+hf'(x)+\frac{h^2}{2}f''(x)+\frac{h^3}{6}f'''(x)+\cdots \]

可得向前差商误差：

\[ f'(x)-\frac{f(x+h)-f(x)}{h} =-\frac{h}{2}f''(x+\theta_1h)=O(h) \]

向后差商误差：

\[ f'(x)-\frac{f(x)-f(x-h)}{h} =\frac{h}{2}f''(x-\theta_2h)=O(h) \]

中心差商误差：

\[ f'(x)-\frac{f(x+h)-f(x-h)}{2h} =-\frac{h^2}{12}\left[f'''(x+\theta_1h)+f'''(x-\theta_2h)\right]=O(h^2) \]

其中 \(0<\theta_1,\theta_2<1\)。

中心差商为何更准

向前、向后差商的一阶误差项没有抵消，所以精度为 \(O(h)\)；中心差商左右对称，偶然地消去了 \(O(h)\) 项，因此精度提升为 \(O(h^2)\)。

1.3 舍入误差与步长选择¶

步长 \(h\) 不是越小越好。截断误差随 \(h\) 减小而减小，但浮点舍入误差中含有除以 \(h\) 的放大效应。实际计算中误差通常先下降后上升。

实践经验

若 \(h\) 太大，Taylor 截断误差占主导。
若 \(h\) 太小，\(f(x+h)-f(x)\) 会发生严重相消，舍入误差占主导。
中心差商通常比单边差商更稳定，但同样需要合理步长。

1.4 插值型数值微分¶

设 \(f(x)\) 在节点 \(x_0,x_1,\cdots,x_n\) 处的函数值为 \(y_i=f(x_i)\)，插值多项式为 \(\phi_n(x)\)，可用

\[ f^{(k)}(x)\approx \phi_n^{(k)}(x),\qquad k=1,2,\cdots \]

来近似导数。

插值余项为

\[ R_n(x)=f(x)-\phi_n(x) =\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x) \]

因此微分误差可通过对余项求导估计。特别地，在节点上有

\[ R_n'(x_i)\approx \frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}'(x_i) \]

这给出了插值求导公式的截断误差来源。

2 数值积分问题¶

2.1 求积公式¶

数值积分的目标是计算

\[ I=\int_a^b f(x)\,dx \]

当原函数难以用初等函数表示时，采用有限个节点的函数值构造近似：

\[ I\approx I_n=\sum_{k=0}^n A_k f(x_k) \]

其中 \(x_k\) 是求积节点，\(A_k\) 是求积系数。

2.2 代数精度¶

定义 2.1（代数精度）

若某求积公式对所有次数不超过 \(m\) 的多项式都精确成立，但对某个 \(m+1\) 次多项式不再精确，则称该求积公式具有 \(m\) 次代数精度。

例 2.1：含导数项的求积公式

设

\[ \int_0^1 f(x)\,dx\approx A_0f(0)+A_1f(1)+B_0f'(0) \]

令 \(f(x)=1,x,x^2\) 分别精确成立：

\[ A_0+A_1=1 \]

\[ A_1+B_0=\frac12 \]

\[ A_1=\frac13 \]

得

\[ A_0=\frac23,\qquad A_1=\frac13,\qquad B_0=\frac16 \]

该公式不能精确计算 \(x^3\) 的积分，因此代数精度为 2。

3 插值型求积¶

3.1 基本思想¶

在 \([a,b]\) 上用插值多项式 \(P_n(x)\) 近似 \(f(x)\)：

\[ P_n(x)=\sum_{i=0}^n f(x_i)l_i(x) \]

其中 \(l_i(x)\) 是 Lagrange 基函数。于是

\[ \int_a^b f(x)\,dx \approx \int_a^b P_n(x)\,dx = \sum_{i=0}^n f(x_i)\int_a^b l_i(x)\,dx \]

求积系数为

\[ A_i=\int_a^b l_i(x)\,dx \]

误差为

\[ R(f)=\frac1{(n+1)!}\int_a^b f^{(n+1)}(\xi(x))\omega_{n+1}(x)\,dx \]

3.2 梯形公式¶

取两个端点 \(a,b\) 作一次插值，得到

\[ \int_a^b f(x)\,dx \approx \frac{b-a}{2}[f(a)+f(b)] \]

这称为 梯形求积公式。

3.3 Simpson 公式¶

取 \(a,\dfrac{a+b}{2},b\) 三个节点作二次插值，得到

\[ \int_a^b f(x)\,dx \approx \frac{b-a}{6} \left[ f(a)+4f\left(\frac{a+b}{2}\right)+f(b) \right] \]

这称为 抛物线公式 或 Simpson 公式。

几何直觉

梯形公式用直线替代被积函数；Simpson 公式用抛物线替代被积函数，所以在光滑函数上通常明显更精确。

4 Newton-Cotes 求积公式¶

4.1 等距节点¶

Newton-Cotes 公式使用等距节点：

\[ x_i=a+ih,\qquad i=0,1,\cdots,n,\qquad h=\frac{b-a}{n} \]

令 \(x=a+th,\ t\in[0,n]\)，则求积公式可写为

\[ \int_a^b f(x)\,dx\approx \sum_{i=0}^n A_i f(x_i) \]

其中

\[ A_i=(b-a)C_i^{(n)} \]

\(C_i^{(n)}\) 称为 Newton-Cotes 系数：

\[ C_i^{(n)} =\frac{(-1)^{n-i}}{n\,i!(n-i)!} \int_0^n \frac{t(t-1)\cdots(t-n)}{t-i}\,dt \]

4.2 常见闭型 Newton-Cotes 公式¶

\(n=1\)：梯形公式¶

\[ \int_a^b f(x)\,dx \approx \frac{b-a}{2}[f(a)+f(b)] \]

系数为

\[ \frac12,\quad \frac12 \]

\(n=2\)：Simpson 公式¶

\[ \int_a^b f(x)\,dx \approx \frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right] \]

系数为

\[ \frac16,\quad \frac23,\quad \frac16 \]

\(n=4\)：Boole 公式¶

令 \(h=\dfrac{b-a}{4}\)，节点为 \(a,a+h,a+2h,a+3h,b\)，则

\[ \int_a^b f(x)\,dx \approx \frac{b-a}{90} \left[ 7f(a)+32f(a+h)+12f(a+2h)+32f(a+3h)+7f(b) \right] \]

例 4.1：积分比较

用梯形公式、Simpson 公式和 \(n=4\) Newton-Cotes 公式计算

\[ \int_{0.5}^{1}\sqrt{x}\,dx \]

准确值为

\[ \int_{0.5}^{1}\sqrt{x}\,dx =\frac23x^{3/2}\Big|_{0.5}^{1} \]

课件中给出的相对误差为：

\[ \begin{array}{c|c} \text{方法} & \text{相对误差}\\ \hline \text{梯形公式} & 9.7\times10^{-3}\\ \text{Simpson 公式} & 7.0\times10^{-5}\\ \text{Newton-Cotes }(n=4) & 7.9\times10^{-7} \end{array} \]

5 误差分析¶

5.1 代数精度结论¶

定理 5.1

\(n+1\) 个节点的 Newton-Cotes 求积公式至少具有 \(n\) 次代数精度。

定理 5.2

当 \(n\) 为偶数时，\(n+1\) 个节点的 Newton-Cotes 求积公式至少具有 \(n+1\) 次代数精度。

这解释了 Simpson 公式的特殊性：它有 3 个节点，本应至少 2 次精确；由于 \(n=2\) 为偶数，实际至少具有 3 次代数精度。

5.2 梯形公式误差¶

若 \(f(x)\in C^2[a,b]\)，则

\[ R_1(f) \equiv \int_a^b f(x)\,dx -\frac{b-a}{2}[f(a)+f(b)] =-\frac{(b-a)^3}{12}f''(\eta) \]

其中 \(\eta\in[a,b]\)。

5.3 Simpson 公式误差¶

若 \(f(x)\in C^4[a,b]\)，则

\[ R_2(f) \equiv \int_a^b f(x)\,dx -\frac{b-a}{6}\left[ f(a)+4f\left(\frac{a+b}{2}\right)+f(b) \right] =-\frac{(b-a)^5}{2880}f^{(4)}(\eta) \]

其中 \(\eta\in[a,b]\)。

5.4 一般 Newton-Cotes 误差¶

对 \(n+1\) 个节点的 Newton-Cotes 公式：

若 \(n\) 为偶数且 \(f(x)\in C^{n+2}[a,b]\)，则

\[ R_n(f) = \frac{h^{n+3}}{(n+2)!}f^{(n+2)}(\eta) \int_0^n t^2(t-1)\cdots(t-n)\,dt \]

若 \(n\) 为奇数且 \(f(x)\in C^{n+1}[a,b]\)，则

\[ R_n(f) = \frac{h^{n+2}}{(n+1)!}f^{(n+1)}(\eta) \int_0^n t(t-1)\cdots(t-n)\,dt \]

其中 \(\eta\in[a,b]\)。

高阶 Newton-Cotes 不一定更好

高阶闭型 Newton-Cotes 公式可能出现负权重和数值不稳定。实际计算中常优先使用复合梯形、复合 Simpson 或 Gauss 求积，而不是盲目提高单个 Newton-Cotes 公式的阶数。

7 复化求积方法¶

7.1 复化梯形公式¶

将区间 \([a,b]\) 进行 \(n\) 等分，节点为

\[ x_k=a+kh,\quad k=0,1,\cdots,n,\quad h=\frac{b-a}{n} \]

在每个小区间 \([x_k,x_{k+1}]\) 上采用梯形公式，然后求和，得到 复化梯形公式：

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

7.2 复化 Simpson 公式¶

由于 Simpson 公式需要区间中点，将 \([a,b]\) 进行 \(2n\) 等分，节点间距仍为 \(h\)，但在每两个节点组成的小区间 \([x_{2k-2},x_{2k}]\) 上采用 Simpson 公式，得到 复化 Simpson 公式：

\[ S_n=\frac{h}{3}\left[f(a)+f(b)+4\sum_{k=1}^{n}f(x_{2k-1})+2\sum_{k=1}^{n-1}f(x_{2k})\right] \]

复化公式的优势

复化公式将高阶 Newton-Cotes 的全局风险转化为局部低阶插值的累积，既保证了精度，又避免了高阶公式可能出现的数值不稳定问题。

7.3 误差估计¶

定理 7.1（复化梯形公式误差）

若 \(f(x)\in C^2[a,b]\)，则

\[ R(f,T_n)\equiv\int_a^b f(x)\,dx-T_n=-\frac{b-a}{12}h^2f''(\eta) \]

其中 \(\eta\in(a,b)\)，\(h=(b-a)/n\)。

定理 7.2（复化 Simpson 公式误差）

若 \(f(x)\in C^4[a,b]\)，则

\[ R(f,S_n)\equiv\int_a^b f(x)\,dx-S_n=-\frac{b-a}{2880}h^4f^{(4)}(\eta)=-\frac{b-a}{180}h_1^4f^{(4)}(\eta) \]

其中 \(\eta\in(a,b)\)，\(h=(b-a)/(2n)\)，\(h_1=(b-a)/n\)。

8 逐次分半法¶

8.1 基本思想¶

复化求积需要预先给定合适的步长 \(h\)：\(h\) 太大则精度不够，\(h\) 太小则计算量过大。更重要的是，估计被积函数高阶导数在区间内的上下界往往十分困难。

逐次分半法 通过不断将区间对半细分，并利用前后两次近似值进行 事后误差估计，直到满足精度要求为止。

8.2 梯形公式的逐次分半¶

设已用 \(n\) 等分得到复化梯形值 \(T_n\)，现将每个小区间增加中点，得到 \(T_{2n}\)。两者之间有关系：

\[ T_{2n}=\frac{1}{2}T_n+\frac{h}{2}\sum_{k=1}^{n}f\left(a+\left(k-\frac{1}{2}\right)h\right) \]

其中 \(h=(b-a)/n\) 为原来的步长。该公式表明：计算 \(T_{2n}\) 时，只需对新增加的中间节点求函数值，原有节点结果可完全复用。

8.3 Simpson 公式的逐次分半¶

类似地，对复化 Simpson 公式有：

\[ S_{2n}=\frac{1}{2}S_n+\frac{h}{3}\sum_{k=1}^{2n}f(x_{k-1/2})-\frac{h}{3}\sum_{k=1}^{n}f(x_{2k-1}) \]

8.4 误差的事后估计¶

定理 8.1（梯形事后估计）

复化梯形公式的误差满足

\[ R(f,T_n)\approx c\,h^2 \]

因此细分后

\[ T_{2n}-T_n\approx 3c\left(\frac{h}{2}\right)^2 \]

从而得到事后估计

\[ I-T_{2n}\approx\frac{T_{2n}-T_n}{3} \]

若容许误差为 \(\varepsilon\)，则当 \(|T_{2n}-T_n|<3\varepsilon\) 时可停止计算。

定理 8.2（Simpson 事后估计）

复化 Simpson 公式的误差满足

\[ R(f,S_n)\approx c\,h^4 \]

从而

\[ I-S_{2n}\approx\frac{S_{2n}-S_n}{15} \]

若容许误差为 \(\varepsilon\)，则当 \(|S_{2n}-S_n|<15\varepsilon\) 时可停止计算。

9 Richardson 外推法¶

Richardson 外推是一种通用的 加速收敛 技术。

假设某个真值 \(F^*\) 的近似值 \(F_1(h)\) 与步长 \(h\) 相关，其截断误差可展开为

\[ F^*-F_1(h)=a_1h^{p_1}+a_2h^{p_2}+\cdots+a_kh^{p_k}+\cdots \]

其中 \(0<p_1<p_2<\cdots\)，\(a_i\neq 0\)。

将 \(qh\) 代替 \(h\)，并与原式进行线性组合，消去 \(h^{p_1}\) 项，得到更高阶的近似：

\[ F_2(h)=\frac{F_1(qh)-q^{p_1}F_1(h)}{1-q^{p_1}} \]

此时 \(F_2(h)\) 逼近 \(F^*\) 的误差阶为 \(O(h^{p_2})\)。

依此类推，可定义递推格式：

\[ F_m(h)=\frac{F_{m-1}(qh)-q^{p_{m-1}}F_{m-1}(h)}{1-q^{p_{m-1}}},\quad m=2,3,\cdots \]

每一步外推都将误差阶提高一级。

10 Romberg 求积¶

Romberg 求积是在复化梯形公式基础上，应用 Richardson 外推构造的一种高效算法。

10.1 基本递推¶

复化梯形公式的误差可写成仅含偶次幂的展开：

\[ R(f,T_n)=a_2h^2+a_4h^4+a_6h^6+\cdots \]

记 \(T_0^{(k)}\) 为将区间 \(2^k\) 等分后的复化梯形值（步长 \(h/2^k\)）。取 \(q=1/2\)，\(p_m=2m\)，应用 Richardson 外推：

\[ T_m^{(k)}=\frac{4^mT_{m-1}^{(k+1)}-T_{m-1}^{(k)}}{4^m-1},\quad m=1,2,\cdots;\ k=0,1,2,\cdots \]

其中： - \(T_0^{(k)}\) 为复化梯形值，误差 \(O(h^2)\) - \(T_1^{(k)}\) 实际上等价于复化 Simpson 值，误差 \(O(h^4)\) - \(T_2^{(k)}\) 误差 \(O(h^6)\) - \(T_m^{(k)}\) 误差 \(O(h^{2m+2})\)

10.2 计算表格¶

Romberg 求积的计算过程可排列成如下三角形表格：

\(k\)	\(m=0\) (梯形)	\(m=1\) (Simpson)	\(m=2\) (Cotes)	\(m=3\) (Romberg)
0	\(T_0^{(0)}\)
1	\(T_0^{(1)}\)	\(T_1^{(0)}\)
2	\(T_0^{(2)}\)	\(T_1^{(1)}\)	\(T_2^{(0)}\)
3	\(T_0^{(3)}\)	\(T_1^{(2)}\)	\(T_2^{(1)}\)	\(T_3^{(0)}\)

计算顺序按列进行：先算第 1 列的梯形值，再依次用外推公式计算右侧各列。

Romberg 的实用价值

Romberg 求积仅需计算梯形序列（函数值可复用），通过简单的代数组合即可迅速提高精度，是数值积分中效率最高的经典算法之一。

例 10.1：Romberg 求积实例

用 Romberg 求积计算

\[ I=\int_0^2 x e^x\,dx \]

其解析解为 \(I=e^2+1\approx 8.389056098930649\)。

计算过程如下：

\(k\)	\(m=0\)	\(m=1\)	\(m=2\)	\(m=3\)
0	14.77811220	8.55041317	8.39037296	8.38905873
1	10.10733793	8.40037547	8.38907926
2	8.82711608	8.38978528
3	8.49911798

可见经过 3 次外推后，结果已精确到小数点后 6 位以上。

11 Gauss 型求积公式¶

11.1 基本思想¶

此前介绍的求积公式（Newton-Cotes）均采用 等距节点。一个自然的问题是：若固定节点数目 \(n\)，能否通过 调整节点位置 来获得更高的代数精度？

考虑一般求积公式

\[ \int_a^b f(x)\,dx\approx\sum_{k=1}^{n}A_kf(x_k) \]

其中 \(n\) 个节点位置 \(x_k\) 和 \(n\) 个系数 \(A_k\) 均待定，共 \(2n\) 个未知数。令其对 \(1,x,x^2,\cdots,x^{2n-1}\) 精确成立，可列 \(2n\) 个方程，因此代数精度 至少可达 \(2n-1\)。达到这一精度的求积公式称为 Gauss 型求积公式。

对于一般区间 \([a,b]\)，可通过变换

\[ x=\frac{a+b}{2}+\frac{b-a}{2}t,\quad t\in[-1,1] \]

将其转化到标准区间 \([-1,1]\) 上讨论。

11.2 两个积分点的 Gauss-Legendre 公式¶

在 \([-1,1]\) 上，取 \(n=2\)，要求

\[ \int_{-1}^{1}f(x)\,dx\approx A_1f(x_1)+A_2f(x_2) \]

对任意三次多项式精确成立。

利用正交性条件：选择 \(x_1,x_2\) 使得 \((x-x_1)(x-x_2)\) 与 \(1\) 和 \(x\) 在 \([-1,1]\) 上正交，解得

\[ x_1=-\frac{1}{\sqrt{3}},\quad x_2=\frac{1}{\sqrt{3}} \]

再确定系数得 \(A_1=A_2=1\)。于是

\[ \int_{-1}^{1}f(x)\,dx\approx f\left(-\frac{1}{\sqrt{3}}\right)+f\left(\frac{1}{\sqrt{3}}\right) \]

这就是 两点 Gauss-Legendre 公式，具有 3 次代数精度。

11.3 一般的 Gauss 型求积公式¶

考虑带权积分

\[ I=\int_a^b w(x)f(x)\,dx \]

其中 \(w(x)\geqslant 0\) 为权函数。一般的 Gauss 型求积公式为

\[ \int_a^b w(x)f(x)\,dx\approx\sum_{k=1}^{n}A_kf(x_k) \]

要求对不高于 \(2n-1\) 次的多项式精确成立。

定理 11.1（Gauss 求积的充要条件）

上述求积公式为 Gauss 型的充分必要条件是：以积分节点 \(x_k\ (k=1,2,\cdots,n)\) 为零点的多项式

\[ \omega_n(x)=(x-x_1)(x-x_2)\cdots(x-x_n) \]

与任意次数不超过 \(n-1\) 的多项式 \(q(x)\) 在区间 \([a,b]\) 上关于权函数 \(w(x)\) 正交，即

\[ \int_a^b w(x)q(x)\omega_n(x)\,dx=0,\quad\forall q(x)\in P_{n-1} \]

此时求积系数可通过 Lagrange 插值基函数计算：

\[ A_k=\int_a^b\frac{w(x)\omega_n(x)}{(x-x_k)\omega_n'(x_k)}\,dx \]

定理 11.2（Gauss 求积误差）

设 \(f(x)\in C^{2n}[a,b]\)，则 Gauss 型积分公式的截断误差为

\[ R=\int_a^b w(x)f(x)\,dx-\sum_{k=1}^{n}A_kf(x_k)=\frac{f^{(2n)}(\eta)}{(2n)!}\int_a^b w(x)\omega_n^2(x)\,dx \]

其中 \(a<\eta<b\)。

11.4 常见的 Gauss 求积公式¶

Gauss-Legendre 求积¶

权函数 \(w(x)=1\)，积分区间 \([-1,1]\)。节点为 Legendre 多项式 \(P_n(x)\) 的零点，系数为

\[ A_k=\frac{2(1-x_k^2)}{(nP_{n-1}(x_k))^2} \]

常见节点与系数：

\(n\)	\(x_k\)	\(A_k\)
1	\(0\)	\(2\)
2	\(\pm 0.5773502692\)	\(1\)
3	\(0,\ \pm 0.7745966692\)	\(0.8888888889,\ 0.5555555556\)
4	\(\pm 0.3399810436,\ \pm 0.8611363116\)	\(0.6521451549,\ 0.3478548451\)

Gauss-Laguerre 求积¶

权函数 \(w(x)=e^{-x}\)，积分区间 \([0,+\infty)\)。节点为 Laguerre 多项式 \(L_n(x)\) 的零点。

Gauss-Hermite 求积¶

权函数 \(w(x)=e^{-x^2}\)，积分区间 \((-\infty,+\infty)\)。节点为 Hermite 多项式 \(H_n(x)\) 的零点。

Gauss-Chebyshev 求积¶

权函数 \(w(x)=1/\sqrt{1-x^2}\)，积分区间 \([-1,1]\)。节点为 Chebyshev 多项式 \(T_n(x)\) 的零点：

\[ x_k=\cos\frac{(2k-1)\pi}{2n},\quad k=1,2,\cdots,n \]

此时所有系数相等：\(A_k=\pi/n\)。

12 数值积分专题¶

12.1 振荡函数的积分¶

工程实际中常需计算形如

\[ \int_a^b f(x)\cos\omega x\,dx,\quad \int_a^b f(x)\sin\omega x\,dx \]

的积分，其中 \(f(x)\) 为非振荡函数，而 \(\omega\) 是较大的正数。直接采用普通数值积分方法需要极多的节点。

利用分段线性插值 \(\phi(x)\) 近似 \(f(x)\)，则每个小区间 \([x_i,x_{i+1}]\) 上的积分

\[ \int_{x_i}^{x_{i+1}}\phi(x)\cos\omega x\,dx \]

可通过分部积分得到解析表达式，从而构造专门针对振荡函数的数值积分公式。其误差满足

\[ \left|\int_a^b f(x)\cos\omega x\,dx-\int_a^b\phi(x)\cos\omega x\,dx\right|\leqslant\frac{b-a}{8}h^2\max_{a\leqslant x\leqslant b}|f''(x)| \]

12.2 奇异函数的积分¶

若被积函数在 \([a,b]\) 内有奇异点，例如

\[ \lim_{x\to a^+}f(x)=\infty \]

且 \(f(x)\) 可写成

\[ f(x)=\frac{\phi(x)}{(x-a)^\mu},\quad 0<\mu<1,\quad |\phi(x)|\leqslant M \]

当 \(\mu\geqslant 1\) 时积分本身发散。

处理方法一：域分解 + Taylor 展开

将积分分解为

\[ I=\int_a^{a+\varepsilon}\frac{\phi(x)}{(x-a)^\mu}\,dx+\int_{a+\varepsilon}^b\frac{\phi(x)}{(x-a)^\mu}\,dx \]

对 \(\phi(x)\) 在 \(x=a\) 处进行 Taylor 展开，第一项可解析计算，余项可用普通数值积分方法。

处理方法二：减去奇异性

令

\[ \phi_p(x)=\sum_{k=0}^{p}\frac{(x-a)^k}{k!}\phi^{(k)}(a) \]

则

\[ I=\sum_{k=0}^{p}\frac{(b-a)^{k+1-\mu}}{k!(k+1-\mu)}\phi^{(k)}(a)+\int_a^b\frac{\phi(x)-\phi_p(x)}{(x-a)^\mu}\,dx \]

右端第二项的被积函数在 \([a,b]\) 内无奇异性。

12.3 无限积分域函数的积分¶

考虑

\[ I=\int_a^{\infty}f(x)\,dx \]

假设 \(I\) 为有限值，一个充分条件是存在 \(\mu>0\) 使得

\[ \lim_{x\to+\infty}x^{1+\mu}f(x)=0 \]

截断法：选择充分大的 \(b\)，使得

\[ I=\int_a^b f(x)\,dx+\int_b^{\infty}f(x)\,dx \]

两部分误差均在可控范围内。

变量替换法：令 \(x=1/t\)，则

\[ I=\int_0^{1/a}f\left(\frac{1}{t}\right)\frac{1}{t^2}\,dt \]

转化为有限区间上的积分，可按奇异函数积分方法处理。

Gauss 型公式：Gauss-Laguerre 和 Gauss-Hermite 求积公式可直接处理半无限和无限区间上的积分。

12.4 自适应数值积分¶

自适应积分（Adaptive Quadrature）根据被积函数的局部变化 自动调整步长或子区间划分：

在函数变化剧烈的地方细分小区间（用更密的节点）
在函数平坦的地方使用较宽区间

从而在给定误差容限的前提下 最小化计算量。

13 多重积分¶

13.1 矩形区域上的累次积分¶

对于二重积分，若积分区域为矩形 \(S=\{(x,y)\mid a\leqslant x\leqslant b,\ c\leqslant y\leqslant d\}\)，则可写成累次积分

\[ \iint_S f(x,y)\,dA=\int_a^b\left[\int_c^d f(x,y)\,dy\right]dx \]

于是前述各种数值积分方法可分别应用于内层和外层积分。

13.2 复合 Simpson 公式¶

对矩形区域，先对 \(y\) 方向、再对 \(x\) 方向应用 Simpson 公式，得到

\[ \iint_S f(x,y)\,dA\approx\frac{(b-a)(d-c)}{36}\Big[ f(a,c)+4f\left(\frac{a+b}{2},c\right)+f(b,c) \]

\[ +4\left(f\left(a,\frac{c+d}{2}\right)+4f\left(\frac{a+b}{2},\frac{c+d}{2}\right)+f\left(b,\frac{c+d}{2}\right)\right) \]

\[ +f(a,d)+4f\left(\frac{a+b}{2},d\right)+f(b,d)\Big] \]

更高维的重积分也可类似地用累次积分方法处理。

14 小结¶

数值微分和数值积分都可以由插值思想统一理解：先用简单函数近似 \(f(x)\)，再对近似函数求导或积分。

数值微分 对数据噪声和舍入误差敏感，步长选择很重要。
数值积分 通常比数值微分稳定，因为积分有平滑效应。
梯形公式 简单稳健，Simpson 公式 在光滑函数上精度高。
Newton-Cotes 公式体现了插值求积的一般框架，但高阶公式需要谨慎使用。
复化求积 通过分段低阶插值兼顾精度与稳定性；逐次分半法 配合事后估计可自适应地控制误差。
Richardson 外推 是一种通用加速收敛技术；Romberg 求积 将其与复化梯形结合，效率极高。
Gauss 型求积 通过优化节点位置达到最高代数精度 \(2n-1\)，是节点数受限时的最优选择。
针对 振荡函数、奇异函数、无限区间 等特殊积分问题，需要采用专门的技巧或变换。
多重积分 可归结为累次积分，逐维应用一维求积方法即可。