数值微分与积分

2026-07-22

用离散点的函数值近似连续的微分/积分，通过多项式插值作为桥梁——先从离散点构造插值多项式，再对多项式求导/积分获得数值微分/积分公式。

📑 目录

• 1. 插值方法概览
• 2. 拉格朗日插值
• 3. 厄密插值
• 4. 龙格现象与分段插值
• 5. 三次样条插值
• 6. Akima 插值
• 7. 多维插值
• 8. 分数多项式插值与 Padé 近似
• 9. 数值微分
• 10. 牛顿-柯特斯积分
• 11. 复合积分与变步长积分
• 12. 理查森外推与龙贝格算法
• 13. 高斯积分
• 14. WKB 近似
• 15. 总结对比表

1. 插值方法概览

                        ┌─────────────────────────┐
                        │      插值方法体系        │
                        └───────────┬─────────────┘
                                    │
          ┌─────────────┬───────────┼───────────┬─────────────┐
          │             │           │           │             │
    全局多项式插值   分段插值    有理插值    多维插值    函数逼近
    ┌───┴───┐    ┌───┴───┐    ┌──┴──┐    ┌──┴──┐    ┌──┴──┐
    │拉格朗日│    │三次样条│    │分数 │    │双线性│    │Taylor│
    │ Hermite│    │ Akima │    │Padé │    │Kriging│   │Cheby │
    └───────┘    └───────┘    └─────┘    └──────┘    └──────┘

插值的本质：已知 N+1 个离散点 (x₀, y₀), (x₁, y₁), ..., (x_N, y_N)，构造一个"穿过"所有点的光滑函数 P(x)，用 P(x) 来估计任意位置的函数值。

2. 拉格朗日插值

2.1 核心思想

人话：给 N+1 个点 $(x₀,y₀)...(x_N,y_N)$，找一个 N 次多项式 P(x) 恰好穿过这些点。做法是构造一组"开关函数" ℓ_k(x)，每个 ℓ_k 在自己那个点等于 1，在其他点等于 0。

2.2 插值基函数

第 k 个基函数 $ℓ_k(x)$ 的构造（N+1 个点，k = 0, 1, ..., N）：

$$ ℓ_k(x) = \frac{(x - x₀)}{(x_k - x₀)} · \frac{(x - x_1)}{(x_k - x_1)} · ... · \frac{(x - x_{k-1})}{(x_k - x_{k-1})} · (\frac{(x - x_{k+1})}{(x_k - x_{k+1})} · ... ·\frac{(x - x_{N})}{(x_k - x_{N})} $$

中文解释：分子是 (x - 除了x_k以外的所有x_i) 连乘，分母是 (x_k - 除了x_k以外的所有x_i) 连乘。当 x=x_k 时分子分母完全相同，结果为1；当 x=x_i（i≠k）时，分子中有一个因子 (x_i - x_i)=0，结果为0。这就是"开关"效果。

2.3 插值多项式

拉格朗日插值多项式 = 各点函数值 × 对应基函数的加权和：

$$ P(x) = y₀·ℓ₀(x) + y₁·ℓ₁(x) + ... + y_N·ℓ_N(x) $$

2.4 具体数值例子

两点插值（直线）：已知 (x₀=1, y₀=2), (x₁=3, y₁=6)

ℓ₀(x) = (x - 3)/(1 - 3) = -(x-3)/2
ℓ₁(x) = (x - 1)/(3 - 1) =  (x-1)/2

P(x) = 2·ℓ₀(x) + 6·ℓ₁(x)
     = 2·(3-x)/2 + 6·(x-1)/2
     = (3-x) + 3(x-1)
     = 3 - x + 3x - 3
     = 2x                          ← 还原为直线 y=2x ✓

三点插值（抛物线）：已知 (x₀=-1, y₀=0), (x₁=0, y₁=1), (x₂=1, y₂=0)

ℓ₀(x) = (x-0)(x-1)/((-1-0)(-1-1)) = x(x-1)/2
ℓ₁(x) = (x+1)(x-1)/((0+1)(0-1))   = -(x+1)(x-1) = 1-x²
ℓ₂(x) = (x+1)(x-0)/((1+1)(1-0))   = x(x+1)/2

P(x) = 0·ℓ₀(x) + 1·ℓ₁(x) + 0·ℓ₂(x)
     = 1 - x²                      ← 还原为抛物线 y=1-x² ✓

💡 物理直觉：拉格朗日插值就像用 N+1 个"开关"拼出一个波形。每个基函数是一个"尖峰"——在自己的点等于1、在别人的点等于0。叠加起来后恰好能穿过所有数据点。

2.5 插值余项（误差）

人话：用 P(x) 代替真实函数 f(x) 的误差有多大？

余项公式（ξ 在 x₀ 和 x_N 之间）：

$$ R(x) = f(x) - P(x) = \frac{f^{N+1}(ξ)}{(N+1)!} · (x-x₀)(x-x₁)...(x-x_N) $$

数值理解：

误差 = 第(N+1)阶导数的值 / (N+1)的阶乘 × 所有距离的乘积
       ↑ 表示函数的弯曲程度      ↑ 表示点之间的距离效应

• 函数越"弯"（高阶导数大）→ 误差越大
• 点越多、阶乘越大 → 分母越大 → 误差可能减小
• 但！龙格现象会破坏这个直觉（见第4节）

3. 厄密插值

3.1 核心思想

人话：拉格朗日只用了函数值 yₖ。如果还知道每一点的导数 y'ₖ，就能做得更精确——用 2N+1 次多项式的精度。

3.2 插值条件

已知 N+1 个点，每个点有 两个信息：

• 函数值： $f(x_{k}) = y_{l}$
• 一阶导数： $f'(x_{k}) = y^{'}_{k}$

总共 2N+2 个条件 → 可以确定 2N+1 次多项式！

3.3 Hermite 基函数

构造两类基函数：

• $hₖ(x)$：在第 k 个点 hₖ(xₖ)=1，且 hₖ'(xₖ)=0；在其他所有点 hₖ(xⱼ)=0, hₖ'(xⱼ)=0
• $ĥₖ(x)$：在第 k 个点 ĥₖ(xₖ)=0，且 ĥₖ'(xₖ)=1；在其他所有点 ĥₖ(xⱼ)=0, ĥₖ'(xⱼ)=0

Hermite 插值多项式：

$$ H(x) = \sum y_{k} h_{k}(x) + \sum y_{k}^{'} \hat{h_{k}} $$

↑ 函数值贡献 ↑ 导数值贡献

设拉格朗日插值基函数为 Li(x)（它在 xi处为1，在其他节点处为0）。我们定义两组新的基函数：

1. A_i(x)：负责控制函数值。它保证在 xi处值为1，导数为0；在其他节点处值为0，导数为0。
2. B_i(x)：负责控制导数值。它保证在 xi处值为0，导数为1；在其他节点处值为0，导数为0。

公式如下（这是你需要记住的黄金公式）：

$$ H(x) = \sum_{i = 0}^{n} [y_i \cdot A_{i}(x) + y_{i}^{'} \cdot B_{i}(x)] $$

其中：

$$ A_i (x) = [1 - 2(x-x_i)L_{i}^{'}(x_i)] \cdot [L_{i}(x)]^{2} $$

$$ B_{i}(x) = (x- x_i) \cdot [L_i(x)]^{2} $$

直观理解： 平方项 [Li(x)]^{2} 保证了在其他节点处函数值和导数值都为0（因为平方让一阶导在节点处也为0）。而前面的线性项 1−2... 则是为了“修正”在 xi 处的导数不为0的缺陷，强行将导数调整为0。

3.4 数值例子：两点 Hermite

已知两个端点 x₀, x₁ 的函数值和导数值：

• f(x₀)=y₀, f'(x₀)=m₀
• f(x₁)=y₁, f'(x₁)=m₁

间距: h = x₁ - x₀

H(x) = y₀·φ₀(t) + y₁·φ₁(t) + m₀·h·ψ₀(t) + m₁·h·ψ₁(t)

其中 t = (x - x₀)/h ∈ [0,1]

φ₀(t) = 2t³ - 3t² + 1    (t=0时为1, t=1时为0, 导数均为0)
φ₁(t) = -2t³ + 3t²        (t=1时为1, t=0时为0, 导数均为0)
ψ₀(t) = t³ - 2t² + t      (t=0时导数为1, t=1时导数为0, 值均为0)
ψ₁(t) = t³ - t²            (t=1时导数为1, t=0时导数为0, 值均为0)

💡 为什么 Hermite 更精确？ 它利用了"速度"信息，就像已知位置和速度来推轨迹比只知位置更准确。在已知导数的场合（如实验数据附带变化率），Hermite 插值用更少的点达到更高的精度。

即三次插值

具体的推导

两个节点 $x_{0} , x_{1}$

引入 $t = \frac{x-x_{0}}{x_1 - x_0} , h = x_1 - x_0$

借助上述归一化后 Lagrange插值中 两个最基本的 一次基函数 为：

$$ L_0(x) = \frac{x-x_1}{x_0-x_1} = \frac{h(t-1)}{-h} = 1-t \\ L_1(x) = \frac{x-x_0}{x_1-x_0} = t $$ 对于x求导，而非对t求导！！！（回归关于x的等式然后求导即可）

$$ L_0^{'}(x) = - \frac{1}{h} \\ L_1^{'}(x_1) = \frac{1}{h} $$ 带入公式

$$ A_i (x) = [1 - 2(x-x_i)L_{i}^{'}(x_i)] \cdot [L_{i}(x)]^{2} $$

$$ B_{i}(x) = (x- x_i) \cdot [L_i(x)]^{2} $$

可得到

$$ A_1(x) = [1-2(x-x_1)L_1^{'}(x_1)]\cdot [L_1(x)]^{2} = -2t^{3}+3t^2 \\ \space \\ A_0(x) = [1-2(x-x_0)L_0^{'}(x_0)]\cdot [L_0(x)]^{2} = (1-t)^2 (1+2t) = 2t^3-3t^2+1 \\ \space \\ B_0(x) = (x-x_0)\cdot [L_0(x)]^2 = h \cdot (t^3 - t^2 + t) \\ \space \\ b_1(t) = h \cdot (t^3 - t^2) $$ 就得到了

$$ H(t) = y_0 (2t^3 - 3t^2 +1) + y_1 (-2t^3 + 3t^2) + y_0^{'}\cdot h(t^3 - 2t^2 +t) + y_1^{'}\cdot h(t^3-t^2) $$

4. 龙格现象与分段插值

4.1 龙格现象

人话：虽然插值点越多多项式次数越高，但边缘误差反而爆炸式增长！

龙格函数例子：f(x) = 1/(1 + x²) 在 [-5, 5] 上等距取点：

取5个点 (N=4)：边缘误差 ≈ 0.4
取9个点 (N=8)：边缘误差 ≈ 0.6
取13个点 (N=12)：边缘误差 ≈ 2.5   ← 越加越糟！
取17个点 (N=16)：边缘误差 ≈ 20    ← 爆炸！

误差趋势（概念图）：

 N=4:     误差 ···  0.4  ···
 N=8:     误差 ······ 0.6 ······
 N=12:    误差 ············· 2.5 ·············
 N=16:    误差 ····························· 20 ···
                                    ↑ 边缘剧烈振荡

💡 物理直觉：高次多项式很"硬"——要同时穿过很多点就会在两端剧烈摇摆。真实物理量通常很"光滑"，不需要这么硬的拟合。解决方案：分段——每段用低次多项式！

4.2 分段插值思路

不用：一个14次多项式穿过15个点  →  ❌ 边缘振荡严重

改用：把区间切成14段
      每段两个点 → 14条直线  →  ⚠️ 角点不光滑
      每段三个点 → 7条抛物线 →  ⚠️ 连接处导数不连续
      每段四个点+光滑条件 → 三次样条 ✓

5. 三次样条插值

5.1 核心思想

人话：把整个区间切成 N 段，每段用一个三次多项式。除了要求穿过数据点，还要求相邻段在连接处函数值、一阶导数、二阶导数都相等——这样整条曲线就像一条柔韧的细木条（旧时绘图用的"样条"）自然弯曲。

5.2 构造框架

假设 N+1 个数据点 (x₀,y₀), ..., (x_N,y_N)，共 N 个小区间。

第 j 段（区间 [xⱼ, xⱼ₊₁]，j=0,...,N-1）：

Sⱼ(x) = aⱼ + bⱼ(x-xⱼ) + cⱼ(x-xⱼ)² + dⱼ(x-xⱼ)³

每段 4 个未知系数，N 段共 4N 个未知数。

5.3 条件清单

5.4 具体推导过程

对于一个均匀间距 h 的例子（3个点：x₀=1, x₁=2, x₂=3；y₀=0, y₁=1, y₂=0）：

步骤1：利用连续性条件消去 dⱼ

   S''ⱼ(xⱼ) = 2cⱼ
   S''ⱼ(xⱼ₊₁) = 2cⱼ + 6dⱼ·h = 2cⱼ₊₁ （二阶导数连续）

→ dⱼ = (cⱼ₊₁ - cⱼ) / (3h)    (j=0,...,N-2)

步骤2：利用一阶导数连续性消去 bⱼ

   S'ⱼ(xⱼ₊₁) = bⱼ + 2cⱼ·h + 3dⱼ·h² = bⱼ₊₁ （一阶导数连续）

→ bⱼ = (yⱼ₊₁ - yⱼ)/h - h·(2cⱼ + cⱼ₊₁)/3

步骤3：得到关于 cⱼ 的三对角方程组

   h·cⱼ₋₁ + 4h·cⱼ + h·cⱼ₊₁ = 3[(yⱼ₊₁ - yⱼ)/h - (yⱼ - yⱼ₋₁)/h]
   
   即 cⱼ₋₁ + 4cⱼ + cⱼ₊₁ = (3/h)·[(yⱼ₊₁-yⱼ)/h - (yⱼ-yⱼ₋₁)/h]

5.5 三对角矩阵形式

对于 N+1 个点（N 个区间），关于 c₁, c₂, ..., c_{N-1} 的方程组：

┌         ┐ ┌    ┐   ┌                                  ┐
│ 4 1 0 0 │ │ c₁ │   │ 3[(y₂-y₁)/h - (y₁-y₀)/h]/h      │
│ 1 4 1 0 │ │ c₂ │   │ 3[(y₃-y₂)/h - (y₂-y₁)/h]/h      │
│ 0 1 4 1 │ │ c₃ │ = │      ...                         │
│ 0 0 1 4 │ │ c₄ │   │ 3[(y₅-y₄)/h - (y₄-y₃)/h]/h      │
└         ┘ └    ┘   └                                  ┘

即：  三对角矩阵 × c向量 = 右端向量（由函数值差分决定）

中文解释：这是一个三对角线性方程组（每行最多3个非零元素），系数矩阵对角线全是 4，上下对角线全是 1，可以用追赶法（Thomas算法）在 O(N) 时间内高效求解。

5.6 边界条件

自然边界条件（最常用）：

S''(x₀) = 2c₀ = 0  →  c₀ = 0
S''(x_N) = 2c_N = 0  →  c_N = 0

这样第一行和最后一行的方程变为：

第一行：4c₁ + c₂ = 右端
最后行：c_{N-2} + 4c_{N-1} = 右端

一次微分边界条件：给定两端的一阶导数值 S'(x₀)=y'₀, S'(x_N)=y'_N：

2c₀ + c₁ = (3/h)·[(y₁-y₀)/h - y'₀]
c_{N-1} + 2c_N = (3/h)·[y'_N - (y_N-y_{N-1})/h]

💡 物理直觉：自然边界条件相当于说"样条两端可以自由弯曲、不受力矩约束"——就像一根细木条，两端不加拧紧的力。这是物理上最自然的状态，也是默认选择。

5.7 算法流程

输入: N+1 个点 (x₀,y₀)...(x_N,y_N)
输出: 每段的系数 aⱼ,bⱼ,cⱼ,dⱼ

步骤1: 计算间距 hⱼ = xⱼ₊₁ - xⱼ
步骤2: 构造三对角方程组 → 求 c₀,...,c_N
步骤3: 回代求 bⱼ 和 dⱼ
步骤4: aⱼ = yⱼ（直接已知）

伪代码:
for j = 0 to N-1:
    h[j] = x[j+1] - x[j]
    
for j = 1 to N-1:   // 构造右端项
    rhs[j] = 3/h * ( (y[j+1]-y[j])/h - (y[j]-y[j-1])/h )

solve_tridiagonal()  → c[0],...,c[N]

for j = 0 to N-1:
    d[j] = (c[j+1] - c[j]) / (3*h)
    b[j] = (y[j+1]-y[j])/h - h*(2*c[j]+c[j+1])/3
    a[j] = y[j]

6. Akima 插值

6.1 核心思想

人话：也是分段三次插值，但不解方程组！直接用局部数据估算导数，然后用 Hermite 两点插值公式拼出曲线。

Akima 插值的优势：

• ✅ 计算极快（不需要求解方程组）
• ✅ 天然局部性——修改一个点只影响相邻几段
• ✅ 有效避免振荡

6.2 一阶导数的估计

对于点 i，用相邻 4 个点的差分斜率做加权平均得到导数估计：

sᵢ = (yᵢ₊₁ - yᵢ) / (xᵢ₊₁ - xᵢ)     ← 右侧斜率
sᵢ₋₁ = (yᵢ - yᵢ₋₁) / (xᵢ - xᵢ₋₁)   ← 左侧斜率

y'ᵢ = (|sᵢ₊₁ - sᵢ|·sᵢ₋₁ + |sᵢ₋₁ - sᵢ₋₂|·sᵢ) / (|sᵢ₊₁ - sᵢ| + |sᵢ₋₁ - sᵢ₋₂|)

权重逻辑：

导数 = 左侧斜率 × 右侧变化剧烈度 + 右侧斜率 × 左侧变化剧烈度
       ────────────────────────────────────────────────────────
                    左右变化剧烈度之和

- 右侧变化剧烈 → 更多信赖左侧信息
- 左侧变化剧烈 → 更多信赖右侧信息
- 左右都平滑 → 平均

6.3 防止除零

当分母 |sᵢ₊₁ - sᵢ| + |sᵢ₋₁ - sᵢ₋₂| = 0 时（连续三个斜率相同）：

y'ᵢ = (sᵢ₋₁ + sᵢ) / 2     ← 直接取平均

6.4 构造插值

得到每个数据点的一阶导数后，每段 [xᵢ, xᵢ₊₁] 用两点 Hermite 插值公式：

已知: f(xᵢ)=yᵢ, f(xᵢ₊₁)=yᵢ₊₁, f'(xᵢ)=y'ᵢ, f'(xᵢ₊₁)=y'ᵢ₊₁

Sᵢ(x) = yᵢ·φ₀(t) + yᵢ₊₁·φ₁(t) + y'ᵢ·h·ψ₀(t) + y'ᵢ₊₁·h·ψ₁(t)
其中 t = (x-xᵢ)/h, h = xᵢ₊₁-xᵢ
      φ₀,φ₁,ψ₀,ψ₁ 为 Hermite 基函数（与 §3.4 相同）

💡 算法思路：Akima = 局部导数加权估计 + Hermite 分段拼接。省掉了样条插值中解三对角方程组那步，因此计算量从 O(N) 进一步降低——但这个 O(N) 本身就很快了，Akima 真正的优势是局部性：修改一个数据点后不需要重算整条曲线。

6.5 边界处理

在端点处需要额外的数据来估算导数：

• 左端点：用最左 3 个点的二次外推估计导数
• 右端点：用最右 3 个点的二次外推估计导数

7. 多维插值

7.1 双线性插值

人话：二维表格数据，先沿 x 方向插值，再沿 y 方向插值（或反过来），结果等价。

场景：已知网格四点 f(x₀,y₀), f(x₁,y₀), f(x₀,y₁), f(x₁,y₁)，求中间点 f(x,y)。

策略1（逐方向插值）：

Step 1: 在 y=y₀ 这行沿 x 插值
        f(x,y₀) = f(x₀,y₀)·(x₁-x)/Δx + f(x₁,y₀)·(x-x₀)/Δx

        在 y=y₁ 这行沿 x 插值
        f(x,y₁) = f(x₀,y₁)·(x₁-x)/Δx + f(x₁,y₁)·(x-x₀)/Δx

Step 2: 用上面两个结果沿 y 插值
        f(x,y) = f(x,y₀)·(y₁-y)/Δy + f(x,y₁)·(y-y₀)/Δy

        Δx = x₁-x₀,  Δy = y₁-y₀

面积权重理解（另一种等价的看法）：

把点 (x,y) 到四个角连成四个小矩形：
f(x,y) = f₁₁·A₁₁ + f₁₀·A₁₀ + f₀₁·A₀₁ + f₀₀·A₀₀

其中 A 是 (x,y) 到对角的矩形面积占整个网格面积的比例
权重之和 = 1（距离归一化到单位1）

7.2 推广到三维

三维插值 = 先做二维插值得到若干平面上的值，再沿第三方向插值。方法可递推。

7.3 散点数据方法

💡 延伸思考：均匀网格的双线性插值在图像缩放中无处不在（cv2.resize 的 INTER_LINEAR 就是它），因为计算极快且效果可接受。

8. 分数多项式插值与 Padé 近似

8.1 有理函数插值

人话：普通多项式遇到函数有极点（分母为零）时会很糟糕。用分子/分母形式的有理函数 R(x) 可以自然描述极点行为——分母的零点就是极点。

形式：

R(x) = P(x)/Q(x)，其中 P 是 u 次多项式，Q 是 v 次多项式

通过 m+1 = u+v+1 个数据点来确定系数。

8.2 Padé 近似

人话：给定函数的 Taylor 展开，找一个分子 u 次 + 分母 v 次的有理函数，使得有理函数展开后的前 u+v 项与 Taylor 展开完全相同。

构造（令 b₀=1）：

函数 f(x) 的 Taylor 级数：
f(x) = a₀ + a₁·x + a₂·x² + a₃·x³ + ...

Padé 近似 R(x) = (p₀ + p₁x + ... + p_u x^u) / (1 + b₁x + ... + b_v x^v)

要求：R(x) 展开后前 u+v+1 项系数与 Taylor 展开一致

系数求解 — 线性方程组：

对于 k ≥ u+1:
a_k + a_{k-1}·b₁ + a_{k-2}·b₂ + ... + a_{k-v}·b_v = 0

这是关于 b₁,b₂,...,b_v 的方程组（v 个方程解 v 个未知数）

解出 b_i 后回代得到 p₀,...,p_u:
p₀ = a₀
p₁ = a₁ + a₀·b₁
p₂ = a₂ + a₁·b₁ + a₀·b₂
...

数值例子：f(x) = eˣ 的 Padé [1,1] 近似

Taylor: eˣ = 1 + x + x²/2 + x³/6 + ...

设 R(x) = (p₀ + p₁x) / (1 + b₁x)

方程（k=2）: a₂ + a₁·b₁ = 0
→ 1/2 + 1·b₁ = 0
→ b₁ = -1/2

回代：p₀ = a₀ = 1
      p₁ = a₁ + a₀·b₁ = 1 + 1·(-1/2) = 1/2

因此 R(x) = (1 + x/2) / (1 - x/2)
         = (2+x)/(2-x)

验证：展开 R(x) = 1 + x + x²/2 + x³/4 + ...
      其中前3项与 eˣ 完全相同！

💡 为什么 Padé 好？ Padé 在极点附近表现远优于 Taylor 截断。例如 1/(1-x) 在 x≈1 附近，Taylor 需要很多项才收敛，而 Padé [0,1] = 1/(1-x) 直接就对了。

9. 数值微分

9.1 核心思想

人话：从拉格朗日插值多项式 P(x) 出发，直接对 P(x) 求导得到数值微分公式。P(x) 是对 f(x) 的近似，所以 P'(x) 就是对 f'(x) 的近似。

9.2 两点公式（前向差分）

从两点插值推导：已知 f(x) 和 f(x+h) 在 $x_0$ 的取值

插值: P(x) = f(x₀)·(x₁-x)/h + f(x₁)·(x-x₀)/h
求导: P'(x) = [-f(x₀) + f(x₁)] / h

f'(x₀) ≈ [f(x₀+h) - f(x₀)] / h     ← 前向差分 (Forward Difference)

误差：O(h) —— 一阶精度

9.3 三点公式

三点中点公式（误差 O(h²)，最常用）：

已知：f(x-h), f(x), f(x+h)

f'(x) ≈ [f(x+h) - f(x-h)] / (2h)       ← 中心差分，O(h²)

三点端点公式（用于边界）：

f'(x₀) ≈ [-3f(x₀) + 4f(x₁) - f(x₂)] / (2h)      ← 左端点，O(h²)
f'(x₂) ≈ [f(x₀) - 4f(x₁) + 3f(x₂)] / (2h)       ← 右端点，O(h²)

数值验证：f(x) = x³，在 x=1 处，h=0.1

真实值: f'(1) = 3

中点公式: [f(1.1) - f(0.9)] / 0.2 = (1.331 - 0.729) / 0.2 = 3.01
                                    误差 ≈ 0.01 = O(h²) ✓

前向差分: [f(1.1) - f(1)] / 0.1 = (1.331 - 1) / 0.1 = 3.31
                                    误差 ≈ 0.31 = O(h) ✓

📊 对比：

- 前向/后向差分 O(h)：用 2 点，精度一般，适合简单估算
- 中心差分 O(h²)：用 3 点，精度高一阶，最常用
- 五点公式 O(h⁴)：用 5 点，精度更高，但需更多数据

9.4 五点公式

五点中点公式：

f'(x₀) ≈ [f(x₀-2h) - 8f(x₀-h) + 8f(x₀+h) - f(x₀+2h)] / (12h)

权重系数：[-1, +8, 0, -8, +1] / (12h)，关于中心点反对称。

五点端点公式：

f'(x₀) ≈ [-25f(x₀) + 48f(x₁) - 36f(x₂) + 16f(x₃) - 3f(x₄)] / (12h)

系数表（对于等距点 x₀, x₀+h, x₀+2h, x₀+3h, x₀+4h）：

9.5 二阶微分

三点公式（从二阶泰勒展开推导）：

f''(x) ≈ [f(x+h) - 2f(x) + f(x-h)] / h²       ← 中心差分，O(h²)

五点二阶微分公式：

f''(x) ≈ [-f(x+2h) + 16f(x+h) - 30f(x) + 16f(x-h) - f(x-2h)] / (12h²)

系数：[-1, +16, -30, +16, -1] / (12h²)，关于中心点对称。

==除了从插值多项式求导，还可以先做样条插值 → 对样条函数直接解析求导，这样得到的导数在整个区间连续且光滑。==

样条插值微分的优雅解决方案：

1. 全局拟合：我们先用上一节学到的三次样条插值（Cubic Spline），将所有数据点拟合成一条全局光滑的曲线。 在每一个区间上，我们都得到了一个极其平滑的解析函数： $S_j(x) = a_j + b_j(x-x_j) + c_j(x-x_j)^2 + d_j(x-x_j)^3$
2. 解析求导：因为我们已经求出了所有的系数 $a_j, b_j, c_j, d_j$，所以我们可以直接在数学上对这个三次多项式求导，得到一阶和二阶导数： $S'_j(x) = b_j + 2c_j(x-x_j) + 3d_j(x-x_j)^2$ $S''_j(x) = 2c_j + 6d_j(x-x_j)$

这样做的巨大优势：

• 全局一阶连续性：在所有数据点的拼接处，样条插值本身就强制要求了一阶和二阶导数的连续性。因此，你求出来的导数曲线 $S'(x)$ 将是一条完全无缝、平滑的连续曲线，而不会像普通差分法那样在点与点之间出现不连续的跳跃。
• 天然的去噪抗噪性：样条插值在求解三对角方程组时，相当于把每一个点的变动通过全局约束“平摊”到了整根铁丝上，它具有物理上最小弯曲能量的限制，因此它对局部的噪声点具有极强的过滤作用。求导出来的结果远比直接做差分要稳定得多！

10. 牛顿-柯特斯积分

10.1 核心思想

人话：从拉格朗日插值多项式出发，直接对整个插值多项式求定积分。

∫f(x)dx ≈ ∫P(x)dx。因为多项式积分有解析公式，所以得到的就是加权和形式的数值积分公式。

数值积分！！！ 工程计算中 很核心的广泛应用的领域

News-Cotes Formulas

10.2 通用形式

对于一个复杂的连续函数 f(x) 无法直接求出原函数 就利用一个好积分的函数去近似 f(x) 然后对于近似的函数求取积分 —— 好的函数 就是 多项式（插值前面的方法）

设我们在等距分点 $x_k = a + k \cdot h$ (其中 步长为 $h = \frac{b-a}{n}$ ) 构造一个 n次的的拉格朗日多项式 $P_{n}(x)$

$$ f(x) \sim P_{n}(x) = \sum_{k = 0}^{n} f(x_k) L_k(x) $$ L(x) 就是Lagrange 基函数：

$$ L_{k}(x) = \prod_{i=0, i \ne k}^{n} \frac{x- x_i}{x_k - x_i} $$ 现在 我们对这个多项式 在 范围 a 到 b上 进行求和积分：

$$ \int_{a}^{b} f(x) dx \approx \int_{a}^{b} P_n(x) dx = \int_{a}^{b} \sum_{k=0}^{n} f(x_k) L_k(x) dx $$

由于积分是线性算子，我们可以将积分号与求和号对调，把已知的节点值 $f(x_k)$ 提出来：

$$ \int_{a}^{b} f(x) dx \approx \sum_{k=0}^{n} f(x_k) \underbrace{\int_{a}^{b} L_k(x) dx}_{A_k} $$

此为 N-C通用形式的由来！

积分值完全等价于 各个节点值 $f(x_{k})$ 的加权和

权重 $A_k = \int_{a}^{b} L_k(x) dx$。

我们对其进行自变量代换，令 $x = a + t \cdot h$（其中 $t \in [0, n]$），则 $dx = h \cdot dt$

$$ A_k = h \cdot \int_{0}^{n} \prod_{i=0, i \neq k}^{n} \frac{t - i}{k - i} dt = (b-a) \cdot \underbrace{\frac{1}{n} \int_{0}^{n} \prod_{i=0, i \neq k}^{n} \frac{t - i}{k - i} dt}_{C_k^{(n)}} $$

这里，$C_k^{(n)}$ 就是著名的柯特斯系数（Cotes Coefficients）。它满足：

• 与 $f(x)$ 毫无关系，只与插值阶数 $n$ 以及节点索引 $k$ 有关。
• 对称性：$C_k^{(n)} = C_{n-k}^{(n)}$（这是由基函数关于中心点对称的代数结构决定的）。
• 归一性：$\sum_{k=0}^{n} C_k^{(n)} = 1$（因为如果 $f(x) = 1$，多项式插值是完全精确的，此时积分值为 $(b-a)$，代入上式必然要求系数和为 1）。

10.3 常用低阶公式

(1) 梯形公式（n=1，2点）

步长 $h = b-a$，节点为 $x_0 = a, x_1 = b$。

拉格朗日基函数为： $L_0(x) = \frac{x - x_1}{x_0 - x_1} = \frac{b - x}{b - a}$ $L_1(x) = \frac{x - x_0}{x_1 - x_0} = \frac{x - a}{b - a}$

求它们的积分： $A_0 = \int_{a}^{b} \frac{b - x}{b - a} dx = \frac{1}{b-a} \left[ bx - \frac{1}{2}x^2 \right]_{a}^{b} = \frac{(b-a)^2}{2(b-a)} = \frac{b-a}{2}$ 同理，由于对称性， $A_1 = \frac{b-a}{2}$。

所以，我们极其自然地得到了梯形公式： $\int_{a}^{b} f(x) dx \approx \frac{b-a}{2} [f(a) + f(b)]$

柯特斯系数为 $[C_0^{(1)}, C_1^{(1)}] = [\frac{1}{2}, \frac{1}{2}]$

区间 [x₀, x₁]，h = x₁-x₀
积分 ≈ h × [f(x₀) + f(x₁)] / 2
系数: C₀^(1) = C₁^(1) = 1/2

几何意义：用梯形面积代替曲线下面积。

误差：$O(h³·f''(ξ))$，有 1 次代数精度（对一次多项式精确）。

具体数值：∫₀¹ x² dx，h=1

梯形: 1 × [0² + 1²]/2 = 0.5
真实值: 1/3 ≈ 0.333
误差 ≈ 0.167

(2) Simpson 公式（n=2，3点，抛物线法则）

步长 $h = \frac{b-a}{2}$，三个节点分别为：$x_0 = a$，$x_1 = \frac{a+b}{2}$，$x_2 = b$。

为了简化积分计算，我们做无损平移，让中点 $x_1$ 落在原点 $0$ 上。 此时，三个节点变为：$-h, 0, h$。 我们来写出 $L_1(x)$ 的表达式：

$L_1(x) = \frac{(x - x_0)(x - x_2)}{(x_1 - x_0)(x_1 - x_2)} = \frac{(x + h)(x - h)}{(0 + h)(0 - h)} = \frac{h^2 - x^2}{h^2}$

对 $L_1(x)$ 进行在 $[-h, h]$ 上的积分：

$A_1 = \int_{-h}^{h} \frac{h^2 - x^2}{h^2} dx = \frac{1}{h^2} \left[ h^2 x - \frac{1}{3}x^3 \right]_{-h}^{h} = \frac{1}{h^2} \left( 2h^3 - \frac{2}{3}h^3 \right) = \frac{4}{3}h$

由于总区间长度为 $b-a = 2h$，所以： $A_1 = \frac{2}{3}(b-a) \implies C_1^{(2)} = \frac{2}{3} = \frac{4}{6}$

由对称性，另外两个系数 $A_0 = A_2$。 又因为所有系数和必须为 $b-a = 2h$，所以： $A_0 + A_2 = 2h - \frac{4}{3}h = \frac{2}{3}h \implies A_0 = A_2 = \frac{1}{3}h = \frac{b-a}{6}$ 即 $C_0^{(2)} = C_2^{(2)} = \frac{1}{6}$。

==Simpson 公式诞生：==

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

区间 [x₀, x₁, x₂]，h = (b-a)/2
积分 ≈ (b-a)/6 × [f(x₀) + 4f(x₁) + f(x₂)]
系数: C₀^(2) = 1/6, C₁^(2) = 2/3, C₂^(2) = 1/6
等价于: h/3 × [1, 4, 1] · [f₀, f₁, f₂]ᵀ

具体数值：

$\int_{0}^{2} x^{2} dx ，h=1$

三个点: x₀=0, x₁=1, x₂=2
f₀=0, f₁=1, f₂=4

Simpson: 2/6 × [0 + 4×1 + 4] = (2/6)×8 = 8/3 ≈ 2.667
真实值: 8/3 ≈ 2.667  ← 精确！(x²是2次多项式，Simpson有3次精度)

误差 ： $O(h⁵·f⁽⁴⁾(ξ))$，有 3 次代数精度。

代数精度（Degree of Precision）：

定义：如果一个数值积分公式对于任意次数不超过 $m$ 的多项式都能完全精确地成立，而对于 $m+1$ 次多项式不精确，则称该公式具有 $m$ 次代数精度。

按照常理，我们用 $n$ 次拉格朗日多项式去近似函数并做积分，得到的公式应该刚好能对 $n$ 次多项式精确，即具有 $n$ 次代数精度。

• 比如梯形公式（$n=1$）对 $f(x)=x$ 精确，代数精度为 $1$。

但是，Simpson公式（$n=2$）创造了一个数学奇迹！ 它使用的是二次插值（抛物线），理论上只能对二次多项式精确。但是，它对任意三次多项式 $f(x) = x^3$ 居然也是完美精确的！ 它的代数精度是 3

为什么会平空多出一阶精度？

我们用泰勒展开在对称中心 $x_1$ 处对误差项进行分析。

由于 Simpson 公式的节点在几何上是完全对称的，当我们把积分误差展开成泰勒级数时，所有奇数阶的误差项会因为正负抵消而彻底消失！

具体来说，对于 $f(x) = x^3$，其三阶导数 $f'''(x)$ 是个常数，其四阶导数 $f^{(4)}(x) = 0$。 而 Simpson 公式的截断误差公式为： $R[f] = -\frac{(b-a)^5}{2880} f^{(4)}(\xi)$ 因为 $f^{(4)}(\xi) = 0$，所以对于任何三次多项式，其误差严格为零。

💡 黄金规律：对于等距节点且对称的 Newton-Cotes 公式，当阶数 $n$ 为偶数时，其代数精度总是为 $n+1$。它会无条件向高处赠送一阶精度！

这就是为什么大家极度钟爱偶数点公式（如 Simpson 公式、五点柯特斯公式）的根本原因。

📝 推导练习提示：把拉格朗日 L₂(x) 带入积分 ∫L₂(x)dx，化简即得 Simpson 公式。

(3) 柯特斯公式（n=4，5点）

区间 [x₀, ..., x₄]，h = (b-a)/4
积分 ≈ (b-a)/90 × [7f(x₀) + 32f(x₁) + 12f(x₂) + 32f(x₃) + 7f(x₄)]
系数: ×(b-a): [7/90, 32/90, 12/90, 32/90, 7/90]
       ≈ [0.078, 0.356, 0.133, 0.356, 0.078]

误差：O(h⁷·f⁽⁶⁾(ξ))，有 5 次代数精度。

10.4 牛顿-柯特斯系数表

💡 注意：n 越大不一定越好！n≥8 时柯特斯系数出现负数，导致舍入误差放大（与龙格现象类似的原因）。实用中都用复合低阶公式。

在插值法那一章，我们学过龙格现象（Runge's Phenomenon）：在等距节点下，高阶多项式会在边界处产生剧烈的、病态的往复震荡 [9.4, 10.5]。

由于 Newton-Cotes 的权重 $A_k$ 本质上就是对拉格朗日基函数 $L_k(x)$ 的积分： 当 $n \ge 8$ 时，由于高阶插值在边界处的急剧震荡，基函数在某些区域的值会变得极大且为负数。这直接导致积分出来的柯特斯系数 $C_k$ 出现了负数 [10.3, 10.5]。

为什么系数出现负数是致命的？

设我们在计算机中进行数值积分，每个点上的函数值都有微小的舍入误差 $\epsilon_k$（由计算机双精度浮点数限制产生）。 我们实际计算出的加权和为： $I_{\text{calc}} = \sum_{k=0}^{n} A_k (f(x_k) + \epsilon_k) = \sum_{k=0}^{n} A_k f(x_k) + \sum_{k=0}^{n} A_k \epsilon_k$

根据误差传播理论，最终积分结果的误差方差正比于： $\sum_{k=0}^{n} A_k^2$

• 当所有 $A_k$ 均为正数时（由于 $\sum A_k = b-a$）： 所有的 $A_k$ 都比较小，平方和 $\sum A_k^2$ 也会非常小。误差被平摊并受到了抑制。
• 当某些 $A_k$ 变为负数时： 由于所有系数的和必须依然等于 $(b-a)$，这会导致其余正的系数必须变得非常大（例如，有些系数变成 $+10$，有些变成 $-9$，它们加起来还是 $1$）。 但在计算平方和时，负号消失了： $(+10)^2 + (-9)^2 = 100 + 81 = 181 \gg 1$ ——> 这意味着，输入数据中哪怕只有 $10^{-16}$ 的极其微小的舍入误差，乘以这些巨大的正负权重后，也会在求和过程中被急剧放大，导致最终的积分结果完全失真！

解决方案：复合低阶公式（Composite Rules）

为了避免高阶公式的灾难，我们在实际中绝不使用单段的高阶 Newton-Cotes 公式。 相反，我们会把大区间 $[a, b]$ 劈成 $M$ 个极小的子区间。在每一个子区间上，我们只使用极其稳定、绝对不会产生负系数的低阶公式（如复合梯形公式或复合 Simpson 公式），最后把所有子区间的积分值累加起来。 这样既保证了计算的极度稳定，又可以通过增加分点来无限逼近真实积分值。

11. 复合积分与变步长积分

11.1 复合积分思想

一大段区间用高阶公式不如切成很多小段、每段用低阶公式，然后把结果加起来。这就是"复合"求积公式。

11.2 复合 Simpson 公式

代数拼装推导

假设我们将大区间 $[a, b]$ 等分为 $n$ 个小段（注意：因为 Simpson 公式本身需要 $3$ 个点才能构成一个子抛物线，所以小区间数 $n$ 必须是偶数）。 步长 $h = \frac{b-a}{n}$。我们把这 $n$ 个区间两两成对地结合起来，一共拼成 $M = n/2$ 个双步长区间：

• 第 1 对：$[x_0, x_2]$，它的三个点是 $x_0, x_1, x_2$。
• 第 2 对：$[x_2, x_4]$，它的三个点是 $x_2, x_3, x_4$。
• 第 3 对：$[x_4, x_6]$，它的三个点是 $x_4, x_5, x_6$。
• $\dots$
• 最后 1 对：$[x_{n-2}, x_n]$，它的三个点是 $x_{n-2}, x_{n-1}, x_n$。

在每一对上，我们单独套用一次单步长为 $h$ 的标准 Simpson 公式：

$\int_{x_0}^{x_2} f(x)dx \approx \frac{h}{3} [1 \cdot f(x_0) + 4 \cdot f(x_1) + 1 \cdot f(x_2)]$

$\int_{x_2}^{x_4} f(x)dx \approx \frac{h}{3} [1 \cdot f(x_2) + 4 \cdot f(x_3) + 1 \cdot f(x_4)]$

$\int_{x_4}^{x_6} f(x)dx \approx \frac{h}{3} [1 \cdot f(x_4) + 4 \cdot f(x_5) + 1 \cdot f(x_6)]$

$\dots$

$\int_{x_{n-2}}^{x_n} f(x)dx \approx \frac{h}{3} [1 \cdot f(x_{n-2}) + 4 \cdot f(x_{n-1}) + 1 \cdot f(x_n)]$

最后，我们把这所有的子区间积分全部加起来。 注意看那些相接的边界点（偶数下标点，如 $x_2, x_4, x_6 \dots$）！

• 节点 $x_1$ 只有第 1 对区间用到，权重为 $4$。
• 节点 $x_2$ 既是第 1 对区间的右端点，又是第 2 对区间的左端点。因此，它的权重被加了两次：$1 + 1 = 2$！
• 同理，所有奇数下标的节点（ $x_1, x_3, x_5 \dots$）都正好落在各自抛物线段的内部，权重永远保持为 $4$；
• 所有偶数下标的节点（ $x_2, x_4, x_6 \dots$），除了最首端的 $x_0$ 和最末端的 $x_n$ 之外，全部是两个抛物线的交界点，所以权重全部变成 $1+1=2$。

因此，复合 Simpson 的系数形式被无缝地拼接出来：

$$ \int_a^b f(x)dx \approx \frac{h}{3} \left[ f_0 + f_n + 4\sum_{\text{奇数 } i} f_i + 2\sum_{\text{偶数 } i \neq 0,n} f_i \right] $$

将 [a,b] 等分成 n 个子区间（n 为偶数），每个子区间 [x_{2k-2}, x_{2k}] 用 Simpson 公式：

步长 h = (b-a)/n

∫_a^b f(x)dx ≈ 
h/3 × [ f₀ + f_n + 4(f₁+f₃+...+f_{n-1}) + 2(f₂+f₄+...+f_{n-2}) ]
      首      尾      奇数下标(×4)            偶数下标(×2，除首尾)

系数模式：1 → 4 → 2 → 4 → 2 → ... → 4 → 2 → 4 → 1

具体数值例子： $\int_{0}^{4} x² dx$，分成 4 段 (n=4)，h=1

分点: x₀=0, x₁=1, x₂=2, x₃=3, x₄=4
f值:     0,    1,    4,    9,   16

复合Simpson:
= 1/3 × [0 + 16 + 4×(1+9) + 2×(4)]
= 1/3 × [16 + 40 + 8]
= 1/3 × 64 = 64/3 ≈ 21.333

真实值: 4³/3 = 64/3 ≈ 21.333  ← 精确！

11.3 误差估计

在数值积分中，误差阶通常被写成 $O(h^p)$

假设我们正在用复合 Simpson 求解一个积分。

当我们把划分的段数 $n$ 加倍（意味着网格步长 $h$ 缩小到了原来的一半 $h_{\text{new}} = \frac{1}{2}h_{\text{old}}$）：

因为复合 Simpson 拥有 $O(h^4)$ 的高精度收敛性，其截断误差中包含一个 $h^4$ 的因子。

$$ \text{Error}_{\text{new}} \approx C \cdot \left(\frac{h}{2}\right)^4 = C \cdot \frac{h^4}{16} = \frac{1}{16} \text{Error}_{\text{old}} $$ 这意味你只需要付出双倍的计算量，误差就会瞬间暴跌至原来的 $\frac{1}{16}$（大约降低了一个数量级还要多）！ 这就是高阶复合积分公式在实际工程中展现出的恐怖威力。

11.4 变步长积分

不知道用多大的步长合适，从大步长开始，不断二分，直到相邻两次计算结果足够接近。

不断二分网络 通过“新结果”与“旧结果” 的差值，实时自我估算 当前的绝对误差！

设原本有 $n$ 段，步长为 $h$。它的复合梯形公式结果为：

$T_n = h \left[ \frac{1}{2}f(a) + f(x_1) + f(x_2) + \dots + f(x_{n-1}) + \frac{1}{2}f(b) \right]$

现在，我们把每一段一分为二，步长变为 $h_{\text{new}} = \frac{h}{2}$，段数变成 $2n$。

新增加的那些中点我们记为 $x_{1/2}, x_{3/2}, \dots, x_{n-1/2}$（共 $n$ 个新增中点）。

现在我们写出全新的 $T_{2n}$（它一共有 $2n+1$ 个点，我们依然把老点和新点在括号里剥离开）： $T_{2n} = h_{\text{new}} \cdot \left\{ \underbrace{\left[ \frac{1}{2}f(a) + f(x_1) + \dots + \frac{1}{2}f(b) \right]}_{\text{所有老点组成的项}} + \underbrace{\sum f(\text{所有新中点})}_{\text{所有新点组成的项}} \right\}$

因为 $h_{\text{new}} = \frac{1}{2} h$，我们把括号前面的 $h_{\text{new}}$ 分配进去：

• 第一项变成了： $\frac{1}{2} h \cdot \left[ \frac{1}{2}f(a) + f(x_1) + \dots + \frac{1}{2}f(b) \right]$ 这正是 $\frac{1}{2}T_n$！
• 第二项保持原样： $h_{\text{new}} \cdot \sum f(\text{所有新中点})$。

于是，我们极其自然地得到了这个递推公式 [11.4]：

$$ T_{2n} = \frac{1}{2} T_n + h_{\text{new}} \sum f\left( 新增中点 \right) $$

变步长梯形算法流程：

步骤1: 初始 n=1
    T₁ = (b-a)·[f(a)+f(b)]/2   ← 一个梯形

步骤2: n 翻倍，利用旧结果
    T₂ₙ = Tₙ/2 + h_new × (新增中点的函数值之和)
    新增中点在: a+h/2, a+3h/2, a+5h/2, ...

步骤3: 检查收敛
    如果 |T₂ₙ - Tₙ| < ε  → 停止
    否则 → 回到步骤2

收敛判定：二分后误差 ≈ 原误差 / 4（对梯形法）

$如果 |T₂ₙ - Tₙ| < ε，则近似积分误差 ≈ |T₂ₙ - Tₙ|/3 < ε/3。$

为什么这里平空冒出了一个除以 3？

这个结论叫做==理查森外推（Richardson Extrapolation）估算==，它的推导极其经典：

我们知道复合梯形公式的真实值 $I$ 与其近似值 $T_n$ 之间的关系可以写成（其中 $C$ 是与步长无关的常数）： $(1) \quad I - T_n = C \cdot h^2 + O(h^4)$

当我们将区间数加倍到 $2n$ 时，新步长变为了 $\frac{h}{2}$。新近似值 $T_{2n}$ 与真实值的关系为： $(2) \quad I - T_{2n} = C \cdot \left(\frac{h}{2}\right)^2 + O(h^4) = \frac{1}{4} C \cdot h^2 + O(h^4)$

我们现在用方程 $(1)$ 减去方程 $(2)$，以此消去未知积分真实值 $I$： $T_{2n} - T_n = \left( C \cdot h^2 \right) - \left( \frac{1}{4} C \cdot h^2 \right) = \frac{3}{4} C \cdot h^2$

由此，我们解出了未知的常量 $C \cdot h^2$： $C \cdot h^2 = \frac{4}{3} (T_{2n} - T_n)$

我们把这个解出来的 $C \cdot h^2$ 重新代回到方程 $(2)$ 中，去看看当前最新计算出来的 $T_{2n}$ 距离真实值 $I$ 到底还差多少：

$\text{当前实际误差} = |I - T_{2n}| \approx \frac{1}{4} C \cdot h^2 = \frac{1}{4} \cdot \left[ \frac{4}{3} (T_{2n} - T_n) \right] = \frac{1}{3} |T_{2n} - T_n|$

证毕！ 这个推导简直妙不可言：我们无法得知积分的真实值 $I$，但我们通过两次近似结果的差值 $|T_{2n} - T_n|$，就能极其精准地估算出当前的实际截断误差恰好就是差值的 $\frac{1}{3}$！ [11.4]

因此，只要我们在程序里监测到 $|T_{2n} - T_n| < \epsilon$（比如 $3 \times 10^{-6}$），我们就能百分之百自信地确信：当前最新的估算值 $T_{2n}$ 的实际误差绝对不会超过 $\frac{\epsilon}{3}$（即 $10^{-6}$）！我们可以立刻安全地终止程序，输出结果。这就是变步长数值积分的精髓所在。

12. 理查森外推与龙贝格算法

12.1 理查森外推思想

人话：一个近似量的误差可以写成 h² 的级数，那么用两个不同步长计算出的结果，可以"消除"最低阶误差项，得到更高精度的近似。

12.2 外推原理

数值微分/积分公式的误差可以写成 h 的偶次幂级数：

D(h) = D_exact + A·h² + B·h⁴ + C·h⁶ + ...
                     ↑ 用两个不同h消掉这项

消去 h² 项：用步长 h 和 h/2 各算一次

D_exact ≈ D(h/2) + [D(h/2) - D(h)] / 3          ← 精度提升两阶！

或者记作: D^(1) = (4·D(h/2) - D(h)) / 3

这就是理查森一阶外推公式！

$$ D^{(1)} = \frac{4 \cdot D\left(\frac{h}{2}\right) - D(h)}{3} = D\left(\frac{h}{2}\right) + \frac{D\left(\frac{h}{2}\right) - D(h)}{3} $$

• 精度跃升：原来的 $D(h)$ 和 $D(\frac{h}{2})$ 都只有二阶精度 $O(h^2)$。
• 但通过这个极其简单的线性组合，我们不费吹灰之力，将最低阶误差项 $h^2$ 彻底消灭，得到的 $D^{(1)}$ 精度瞬间暴涨到了 四阶精度 $O(h^4)$ ！

12.3 数值例子（微分外推）

用中心差分计算 f'(1)，f(x)=sin x：

步长 h=0.4:
  D(0.4) = [sin(1.4)-sin(0.6)]/0.8
         = (0.9854-0.5646)/0.8 = 0.5260

步长 h=0.2:
  D(0.2) = [sin(1.2)-sin(0.8)]/0.4
         = (0.9320-0.7174)/0.4 = 0.5365

外推一次: D_exact ≈ D(0.2) + [D(0.2)-D(0.4)]/3
                 = 0.5365 + 0.0105/3 = 0.5400

真实值: cos(1) = 0.5403    ← 外推后精度大幅提升！
原始误差 0.0142 → 0.0003

12.4 Romberg 积分

理查德外推 可以将精度提升两阶的好用的东西 那么 Romberg积分 就是把这一手段无限嵌套、循环使用的终极策略！！！

理查森外推 + 复合梯形公式 = Romberg 积分。

梯形法逐次二分 → 用不同步长的结果外推 → 精度跃升。

Romberg算法的直观理解：

复合梯形公式 $T_0(h)$ 的精度是 $O(h^2)$

如果我们把 $T_0(h)$ 和 $T_0(\frac{h}{2})$ 进行一次理查森外推，消去 $h^2$ 项，得到的新结果其精度为 $O(h^4)$。在数学上，这个结果恰好等价于复合 Simpson 公式 $T_1(h)$ 的值 ！

同理，既然我们有了好几个不同步长的 Simpson 值 $T_1(h)$ 和 $T_1(\frac{h}{2})$（精度均为 $O(h^4)$），那我们能不能对它们再做一次外推，消去 $h^4$ 项？ 可以！外推后，精度再次暴涨两阶，达到 $O(h^6)$，这恰好等价于复合 Cotes 公式 $T_2(h)$。

以此类推，我们对 Cotes 的值再做外推，就能得到拥有 $O(h^8)$ 恐怖精度的 Romberg 值 $T_3(h)$ ！

Romberg 表格（T-表）：

第0列(T₀): 复合梯形     O(h²)
第1列(T₁): Simpson       O(h⁴)   = (4×T₀(h/2) - T₀(h)) / 3
第2列(T₂): Cotes         O(h⁶)   = (16×T₁(h/2) - T₁(h)) / 15
第3列(T₃): Romberg       O(h⁸)   = (64×T₂(h/2) - T₂(h)) / 63

核心递推公式

$$ T_{j}^{k} = \frac{4^{j} · T_{j-1}^{k+1} - T_{j-1}^{k} }{(4^{j} - 1)} $$ j ：代表列 精度阶数 也就是 外推了多少次

k：代表行 区间划分次数 即二分了多少次

j = 0 j = 1 j = 2 j = 3 梯形列 Simpson列 Cotes列 Romberg列 (原始计算) (一阶外推) (二阶外推) (三阶外推) ─────────────────────────────────────────────────────────────

k = 0 (1段) T_0^(0) ──┐ │ k = 1 (2段) T_0^(1) ──┴───→ T_1^(0) ──┐ │ │ k = 2 (4段) T_0^(2) ──┴───→ T_1^(1) ──┴───→ T_2^(0) ──┐ │ │ │ k = 3 (8段) T_0^(3) ──┴───→ T_1^(2) ──┴───→ T_2^(1) ──┴───→ T_3^(0)

现在你再去看这个外推公式： $T_{\color{blue}j}^{({\color{red}k})} = \frac{4^j \cdot T_{{\color{blue}j-1}}^{({\color{red}k+1})} - T_{{\color{blue}j-1}}^{({\color{red}k})}}{4^j - 1}$

就会发现，要计算一个新位置 $T_j^{(k)}$，它只依赖于它左边那一列（$j-1$）的两个相邻数3：

• 一个是它左下方的数： $T_{j-1}^{(k+1)}$
• 一个是它正左方的数： $T_{j-1}^{(k)}$

我们要算 Simpson 列的第一项 $T_1^{(0)}$ ( $j=1, k=0$)：

我们需要它左下方的 $T_0^{(1)}$ 和正左方的 $T_0^{(0)}$： $T_1^{(0)} = \frac{4 \cdot T_0^{(1)} - T_0^{(0)}}{3}$

我们要算 Cotes 列的第一项 $T_2^{(0)}$ ( $j=2, k=0$) [12.5]：

我们需要它左下方的 $T_1^{(1)}$ 和正左方的 $T_1^{(0)}$： $T_2^{(0)} = \frac{16 \cdot T_1^{(1)} - T_1^{(0)}}{15}$

可以来拆解下不同列 就是 不同 j 时候的具体的代数形式：

1、j = 1 （k = 0 做一次二分）从梯形外推到 Simpson：

带入数值 可以得到

$T_{1}^{(k)} = \frac{4\cdot T_0(h/2)- T_0(h)}{3}$

用来消去 $h^2$ 项的标准理查森外推公式

当 $j=2$ 时（从 Simpson 外推到 Cotes）：

因为此时我们要消去的是 $h^4$ 项，它的误差变化比例是 $\left(\frac{1}{2}\right)^4 = \frac{1}{16}$。

所以我们要用系数 $16$ 来消元 [12.4]：

$T_2^{(k)} = \frac{4^2 \cdot T_1^{(k+1)} - T_1^{(k)}}{4^2 - 1} = \frac{16 \cdot T_1(h/2) - T_1(h)}{15}$

当 $j=3$ 时（从 Cotes 外推到 Romberg）：

消去 $h^6$ 项，误差变化比例是 $\left(\frac{1}{2}\right)^6 = \frac{1}{64}$

$T_3^{(k)} = \frac{4^3 \cdot T_2^{(k+1)} - T_2^{(k)}}{4^3 - 1} = \frac{64 \cdot T_2(h/2) - T_2(h)}{63}$

数学的和谐之美在此处展现得淋漓尽致。所有的外推步骤都被优雅地统一在了一个极其简单的公式 $T_j^{(k)}$ 中。在写代码时，只需要两层循环，就能把整张 Romberg 表源源不断地算出来。

12.5 Romberg 完整数值例子

计算 ∫₀¹ 1/(1+x) dx = ln 2 ≈ 0.693147

步骤1 — 逐次二分梯形:

n=1, h=1:
  T₀^(0) = 1·[1 + 1/2]/2 = 0.750000

n=2, h=0.5:
  T₀^(1) = 0.5/2·[1 + 2×(1/1.5) + 1/2]
         = 0.25·[1 + 1.3333 + 0.5] = 0.708333

n=4, h=0.25:
  T₀^(2) = 0.25/2·[1 + 2×(1/1.25+1/1.5+1/1.75) + 1/2]
         = 0.125·[1 + 4.0762 + 0.5] = 0.697024

n=8, h=0.125:
  T₀^(3) = ... = 0.694122

步骤2 — Romberg 外推:

Simpson列:
  T₁^(0) = (4×0.708333 - 0.75)/3 = 0.694444
  T₁^(1) = (4×0.697024 - 0.708333)/3 = 0.693254
  T₁^(2) = (4×0.694122 - 0.697024)/3 = 0.693155

Cotes列:
  T₂^(0) = (16×0.693254 - 0.694444)/15 = 0.693175
  T₂^(1) = (16×0.693155 - 0.693254)/15 = 0.693148

Romberg列:
  T₃^(0) = (64×0.693148 - 0.693175)/63 = 0.693147 ← 精度极高！

Romberg T-表可视化：

k    梯形 O(h²)    Simpson O(h⁴)   Cotes O(h⁶)    Romberg O(h⁸)
0    0.750000
1    0.708333      0.694444
2    0.697024      0.693254        0.693175
3    0.694122      0.693155        0.693148        0.693147

💡 物理直觉：Romberg 就像"梯度下降"逼近真实值——每一步二分梯形只提供 O(h²) 精度，但两次二分结果组合一下就能跃升到 O(h⁴)，再组合到 O(h⁶)...到最后只用梯形公式算 8 段就能达到 O(h⁸) 精度！

13. 高斯积分

13.1 核心思想

牛顿-柯特斯用 n+1 个等距点达到最多 n+1 次代数精度。但如果允许节点位置自由选择，用同样的 n+1 个点可以达到 2n+1 次代数精度！——这就是高斯积分的精髓。

直观对比：

牛顿-柯特斯 (等距点):  n+1 个点 → n+(n的奇偶修正) 次精度
高斯积分 (最优选点):  n+1 个点 → 2n+1 次精度   ← 翻倍！

核心思想：打破等距的限制！

在前面的 N-C公式中 都有一个默认的限制：节点 x_i 必须是等距离分布的

如果我们固定使用 n+1 个等距节点， 我们最多只能获得 n+1 次代数精度（如果n为偶数，可以利用对称性提升到 n+2 阶精度）

如果我们把这 $n+1$ 个节点 $x_0, x_1, \dots, x_n$ 的位置也当作可以自由调整的未知数，那么：

• 我们有 $n+1$ 个可自由移动的节点位置 $x_i$；
• 我们有 $n+1$ 个对应的求积权重 $A_i$。

一共有 2n+2 自由度 理论上可以写出 2n+2 个代数方程 从而使得 求积公式 对所有次数 不超过 2n+1 的多项式 都达到 绝对精确

13.2 高斯点条件

定理： $x₀,...,x_n$ 是高斯点的充要条件——多项式 $ω(x) = (x-x₀)(x-x₁)...(x-x_n)$ 与所有不超过 n 次的多项式 P_n(x) 在 [a,b] 上正交。

为什么正交能带来高精度？

这里有一个极其优美的代数证明： 对于任意一个最高次数为 $2n+1$ 的多项式 $f(x)$，我们都可以用 $\omega(x)$（它是 $n+1$ 次多项式）去除它。根据多项式带余除法，一定可以写成： $f(x) = q(x)\omega(x) + r(x)$ 其中：

• 商式 $q(x)$ 的次数最高为 $n$ 次。
• 余式 $r(x)$ 的次数最高也为 $n$ 次。

现在，我们对 $f(x)$ 求精确积分： $\int_a^b f(x)dx = \int_a^b q(x)\omega(x)dx + \int_a^b r(x)dx$ 如果 $\omega(x)$ 正交于所有不超过 $n$ 次的多项式（而 $q(x)$ 正好不大于 $n$ 次），那么根据正交性的定义，第一项积分直接归零 [13.2, 13.3]： $\int_a^b q(x)\omega(x)dx = 0$ 因此，真实积分值完全退化为： $\int_a^b f(x)dx = \int_a^b r(x)dx$

接下来，我们用高斯公式来计算 $f(x)$ 的近似积分。

由于高斯节点 $x_0, \dots, x_n$ 正好是 $\omega(x)$ 的零点（即 $\omega(x_i) = 0$），所以： $f(x_i) = q(x_i)\omega(x_i) + r(x_i) = r(x_i)$

因此高斯积分公式给出： $\sum_{i=0}^n A_i f(x_i) = \sum_{i=0}^n A_i r(x_i)$

因为 $r(x)$ 只是一个最高 $n$ 次的多项式，而任何包含 $n+1$ 个节点的插值型求积公式对 $n$ 次多项式都是绝对精确的： $\sum_{i=0}^n A_i r(x_i) = \int_a^b r(x)dx$

把这几步串起来：

$$ 高斯近似值 = \sum_{i=0}^n A_i f(x_i) = \sum_{i=0}^n A_i r(x_i) = \int_a^b r(x)dx = \int_a^b f(x)dx = 真实值 $$ 其中 权重 $A_i$：是这些代表所拥有的“投票权大小”（即话语权比重）

我们要算一段区间的积分（也就是求围成的图形面积），高斯积分就是只挑选 $n+1$ 个特定位置的函数值 $f(x_i)$ 进行加权求和： $\int_a^b f(x) dx \approx A_0 f(x_0) + A_1 f(x_1) + \dots + A_n f(x_n)$

证毕！ 这个证明堪称数学史上最优雅的篇章之一。它表明：只要把节点选在正交多项式的零点上，高阶多项式中高于 $n$ 次的部分就会因为正交性被“自动过滤”掉，使公式精度奇迹般地翻倍！

13.3 高斯-勒让德积分

在 [-1,1] 上，取勒让德多项式 $P_{n+1}(x)$ 的零点作为高斯点。

为什么是勒让德？ 勒让德多项式天然在 [-1,1] 上正交！

勒让德多项式

P₀(x) = 1
P₁(x) = x
P₂(x) = (3x² - 1)/2
P₃(x) = (5x³ - 3x)/2
P₄(x) = (35x⁴ - 30x² + 3)/8

递推: (n+1)·P_{n+1}(x) = (2n+1)·x·P_n(x) - n·P_{n-1}(x)
正交性: ∫^{1}_{-1} P_m(x)·P_n(x) dx = 0   (m ≠ n)
归一化: ∫^{1}_{-1} [P_n(x)]² dx = 2/(2n+1)

求节点：以 n=2（3点）为例

3 个高斯点 = P₃(x) 的 3 个零点：

P₃(x) = (5x³ - 3x)/2 = 0
→ x·(5x² - 3) = 0
→ x₀ = -√(3/5) ≈ -0.774597
  x₁ = 0
  x₂ = √(3/5) ≈ +0.774597

求权重：解矩方程

高斯积分公式： $\int_{-1}^{1} f(x)dx ≈ A₀·f(x₀) + A₁·f(x₁) + A₂·f(x₂)$

要求对 $f(x)=1, x, x², x³, x⁴, x⁵$ 都精确（3 点可达 2×3-1=5 次精度）：

对 f(x)=1:  ∫1 dx = 2 = A₀+A₁+A₂
对 f(x)=x:  ∫x dx = 0 = A₀·(-√0.6) + A₁·0 + A₂·(√0.6)
→ A₀ = A₂

对 f(x)=x²: ∫x²dx = 2/3 = A₀·(0.6) + A₁·0 + A₂·(0.6)
→ 2A₀·0.6 = 2/3,  A₀ = 5/9

→ A₀ = A₂ = 5/9,  A₁ = 2 - 10/9 = 8/9

矩阵形式理解（矩方程）：

┌                     ┐ ┌    ┐   ┌               ┐
│  1      1      1    │ │ A₀ │   │  ∫1·x⁰dx = 2  │
│ x₀      x₁     x₂   │ │ A₁ │   │  ∫1·x¹dx = 0  │
│ x₀²     x₁²   x₂²   │ │ A₂ │ = │  ∫1·x²dx =2/3 │
└                     ┘ └    ┘   └               ┘

代入 x₀=-√0.6, x₁=0, x₂=√0.6，解线性方程组得权重

13.4 高斯-勒让德节点与权重表

13.5 积分区间变换

由于勒让德多项式天然正交于区间 $[-1, 1]$ [13.3]，所有标准的高斯点和权重表都是针对 $[-1, 1]$ 设计的 [13.4]。 如果实际积分区间是 $[a, b]$，我们必须进行一次线性自变量替换。

令 $x = \alpha t + \beta$。我们希望当 $t = -1$ 时 $x = a$，当 $t = 1$ 时 $x = b$。

代入解方程组： $\begin{cases} -\alpha + \beta = a \\ \alpha + \beta = b \end{cases} \implies \alpha = \frac{b-a}{2}, \quad \beta = \frac{a+b}{2}$

由此得到经典的区间映射变换式： $x = \frac{b-a}{2}t + \frac{a+b}{2}$

同时，微分项也需要按比例缩放： $dx = \frac{b-a}{2} dt$ 权重也需要同步缩放 ！

高斯-勒让德节点在 [-1,1] 上，实际积分区间 [a,b] 需要变换：

∫_a^b f(x)dx = (b-a)/2 × ∫_{-1}^{1} f((a+b)/2 + (b-a)t/2) dt

              ≈ (b-a)/2 × Σ Aₖ·f((a+b)/2 + (b-a)xₖ/2)
                                                 k

通俗转化:
变量替换: x = (a+b)/2 + (b-a)t/2
→ dx = (b-a)/2 · dt
→ t∈[-1,1] 时 x∈[a,b]
权重也缩放: Aₖ' = (b-a)/2 × Aₖ

13.6 数值例子（n=2，3点高斯）

计算 ∫₋₁¹ x⁴ dx = 2/5 = 0.4

节点: x₀=-√0.6, x₁=0, x₂=√0.6
权重: A₀=5/9, A₁=8/9, A₂=5/9

高斯积分 = A₀(x₀)⁴ + A₁(0)⁴ + A₂(x₂)⁴
        = 5/9 × 0.36 + 0 + 5/9 × 0.36
        = 10/9 × 0.36 = 0.4  ← 精确！

而Simpson (3个等距点): 
  = 1/3 × [(-1)⁴ + 4×0⁴ + 1⁴] = 2/3 ≈ 0.667  ← 误差巨大！

原因: Simpson对3次以下精确，x⁴是4次 → 有误差
      3点高斯对5次以下都精确！

13.7 带权高斯积分

解决数值积分中的两个超级痛点：

1. 积分区间是无穷大（ $\infty$）怎么办？
2. 函数在端点会爆炸（奇点， $\frac{1}{0}$）怎么办？

公式形式： $∫ ρ(x)·f(x) dx ≈ Σ Aₖ·f(xₖ)$ ，节点 xₖ 是相应正交多项式的零点。

把困难的会爆炸的或者趋于无穷的部分 剥离给 $\rho(x)$ 让可以进行 Guass积分处理的部分 留在f（x）！

💡 延伸思考：高斯-勒让德是"默认选项"、精度极高，但节点不在积分区间的端点（这对某些问题不方便）。高斯-拉盖尔伽利厄米积分天然处理无穷区间，在量子力学中极其重要——例如 ∫₋∞^∞ e^{-x²} f(x) dx 用 Gauss-Hermite 只需几个点就有高精度。

本质上是“正交多项式家族”针对不同物理舞台的量身定制 ！！！

用最好的计算点数 避开所有数学陷阱 算出精确的积分！

14. WKB 近似

14.1 物理背景

WKB（Wentzel-Kramers-Brillouin）方法是量子力学中近似处理势垒穿透（tunneling）问题的半经典方法。一个粒子要穿过势垒的概率可以用一个积分来表示。

14.2 WKB 公式

WKB 方法的物理本质，是将量子波函数 $\psi(x)$ 写成关于普朗克常数 $\hbar$ 的渐近展开式

$$ \psi(x) \sim \exp\left( \frac{i}{\hbar} S(x) \right) $$ 其中 $S(x)$ 对应经典力学中的哈密顿-雅可比作用量。

• 经典极限：当 $\hbar \to 0$ 时，体系的行为完全退化为牛顿经典力学。
• 量子修正：当 $\hbar$ 有限大时，它保留了波的相位信息，从而能够描述经典粒子绝对无法完成的“穿墙术”——量子隧穿

穿透因子（穿透几率）物理意义：代表粒子穿过势垒 $[b, c]$ 的隧穿几率（即成功率）

$$ W = exp(-2 × ∫_b^c \frac{\sqrt{[2M(q)·(V(q)-E₀)]}}{\hbar } dq) $$

周期因子：

$$ T = 2 ∫_a^b \frac{dq}{\sqrt{\frac{2[(E₀-V(q))]}{M(q)}}} $$

寿命：τ = T / W

14.3 公式含义

          V(q)
          ↑
     E₀ ──├─────┐      ← ┌── 势垒（V > E₀，经典禁戒区）
          │     │         │
          │ a   │ b    c  │
    ──────┘     └─────────└──→ q

   a, b: 经典允许区边界（V(a)=V(b)=E₀）
   b, c: 势垒区（V(q) > E₀，粒子需要隧穿）
   [a,b]: 周期 T 的积分区间 → 粒子在阱中的往返时间
   [b,c]: 穿透因子 W 的积分区间 → 决定穿透概率的指数

14.4 数值实现思路

给定表格形式的 V(q) 和 M(q) 离散数据（形变 q 从 0 到 2.5），计算隧穿寿命 τ：

步骤1: 插值 V(q) 和 M(q)（用三次样条获得光滑曲线）
步骤2: 确定临界点 a, b, c（解方程 V(q)=E₀）
        → 用二分法或牛顿法在插值曲线上找根
步骤3: 用高斯积分计算 ∫_b^c √(M(q)(V(q)-E₀)) dq
步骤4: 计算 W = exp(-2×积分结果)
步骤5: 用高斯积分计算 ∫_a^b dq/√((E₀-V(q))/M(q))
步骤6: T = ℏ × 步骤5的结果
步骤7: τ = T / W

数值细节：

• 临界点 a, b, c 通过插值后的 V(q) 曲线用二分法或牛顿法找到
• 在 b 和 c 附近 V(q)≈E₀，被积函数趋于 0 或无穷 → 需小心处理端点奇性
• 被积函数剧烈变化 → 高斯积分比等距公式更高效

💡 物理直觉：W 是量子力学中的穿透因子——势垒越宽越高，W 越小，穿透越不可能。τ 是经典周期除以穿透概率——只有在势垒内来回反射很多次才会偶尔穿透一次，所以寿命 τ 可以比经典周期长很多数量级。这就是为什么某些放射性核素的半衰期可以长达数十亿年。

15. 总结对比表

15.1 插值方法对比

15.2 数值微分公式对比

15.3 数值积分方法对比

15.4 方法选择决策树

需要做什么？
├── 只有离散数据点，需要中间值
│   ├── 点数少（=10）→ 拉格朗日/ Hermite 插值
│   ├── 点数多，要光滑 → 三次样条
│   └── 点数多，要效率 → Akima 插值
│
├── 需要一阶导数
│   ├── 有解析函数 → 中心差分 O(h²)
│   ├── 有离散数据 → 先样条插值再求导
│   └── 边界处 → 三点端点公式
│
├── 需要定积分
│   ├── 被积函数光滑，任意点可求值
│   │   ├── 要绝对高精度 → 高斯积分（n 点有 2n+1 次精度）
│   │   ├── 渐进式提高精度 → Romberg 积分
│   │   └── 简单快速 → 复合 Simpson
│   ├── 只有等距表格数据 → 复合 Simpson/柯特斯
│   ├── 积分区间无限
│   │   ├── [0,∞) 带 e^{-x} 衰减 → Gauss-Laguerre
│   │   └── (-∞,∞) 带 e^{-x²} 衰减 → Gauss-Hermite
│   └── 有奇点或剧烈振荡 → 自适应积分 + 分段
│
└── 需要WKB隧穿计算
    └── 样条插值 + 高斯积分 + 二分法找临界点

15.5 核心概念速记

Akurio 补充整理 笔记结束 · 2026-06-30