5.2 5.2 插值、数值积分与方程求根

5.2 插值、数值积分与方程求根

本节摘要：插值、求积、求根是数值方法的三大常备药。本节讲拉格朗日插值与三次样条的分界线（龙格现象）、复合梯形与辛普森公式的误差阶、自适应积分的机理，以及二分法与牛顿法这对"稳与快"搭档在换热器设计工单中的配合使用。

学习目标

阅读完本节，你应当能够：

1. 解释龙格现象并据此选择多项式插值或样条插值；
2. 用误差阶比较梯形与辛普森公式，说明自适应积分何时启动；
3. 组合二分法与牛顿法实现"快而不断线"的求根器。

插值：过点容易，像样难

热工实验给了十一个温度-黏度标定点，要求一条"能查任意温度"的曲线。最直接的答案是十次多项式插值——精确穿过每个点。问题在点之间：等距节点上的高次多项式会在区间两端剧烈摆振，这就是龙格现象。节点越密，摆振越凶，与直觉完全相反。

import numpy as np

runge = lambda x: 1.0 / (1 + 25*x**2)

def lagrange(nodes, values, xs):
    """拉格朗日插值：基函数加权求和"""
    total = np.zeros_like(xs)
    for i, (xi, yi) in enumerate(zip(nodes, values)):
        basis = np.ones_like(xs)
        for j, xj in enumerate(nodes):
            if i != j:
                basis *= (xs - xj) / (xi - xj)
        total += yi * basis
    return total

xs = np.linspace(-1, 1, 400)
for n in [6, 10, 16]:
    nodes = np.linspace(-1, 1, n+1)
    vals = runge(nodes)
    err = np.max(np.abs(lagrange(nodes, vals, xs) - runge(xs)))
    print(f"等距 {n+1} 点高次插值: 最大误差 {err:.3f}")

# 对照：分段三次样条
from scipy.interpolate import CubicSpline
nodes = np.linspace(-1, 1, 11)
cs = CubicSpline(nodes, runge(nodes))
err_spline = np.max(np.abs(cs(xs) - runge(xs)))
print(f"11 点三次样条:      最大误差 {err_spline:.5f}")
​

等距点数从七加到十七，高次插值误差反而恶化几个数量级；同样十一个点的三次样条误差安静地停在 10 的负 4 次方量级。选型规则：节点少于六七个用多项式没问题；更多点用分段低次（样条），或者干脆用第 2 章的最小二乘——如果数据带噪声，"穿过每个点"本身就是错误目标，插值会把噪声也穿起来。

数值积分：误差阶与自适应

定积分的数值近似有两大门派：单项式区间的牛顿-柯特斯公式（梯形、辛普森）与基于正交多项式的高斯求积。复合梯形的误差与步长平方成正比，复合辛普森与步长四次方成正比——阶数就是"步长减半、误差降多少倍"。

import numpy as np

f = lambda x: np.exp(-x**2)          # 被积函数：误差函数核心
a, b = 0.0, 2.0
ref = 0.882081390762421              # 高精度参考值

def trapezoid(n):
    xs = np.linspace(a, b, n+1)
    h = (b - a) / n
    return h * (f(xs[0])/2 + f(xs[1:-1]).sum() + f(xs[-1])/2)

def simpson(n):
    xs = np.linspace(a, b, n+1)
    h = (b - a) / n
    w = np.ones(n+1); w[1:-1:2] = 4; w[2:-1:2] = 2
    return h/3 * (w * f(xs)).sum()

for n in [10, 20, 40]:
    et = abs(trapezoid(n) - ref)
    es = abs(simpson(n) - ref)
    print(f"n={n:>3}: 梯形误差 {et:.2e}  辛普森误差 {es:.2e}")
​

步长减半时梯形误差约降四倍（二阶）、辛普森约降十六倍（四阶）——输出的比值就是阶数的现场演示。生产中不手动调 n，而是用自适应积分：先在整个区间算一次，再对半拆开算一次，两值之差若小于容差就收工，否则把差的那半段继续拆。工作量自动集中被积函数"陡峭"的区域，平滑处一笔带过。科学计算库的积分函数内部就是这个循环，调用时给容差即可。

另一条经验法则：周期函数在整个周期上积分，梯形公式意外地超高精度（误差按指数衰减）——傅里叶系数计算全靠这条性质。被积函数的特性比公式的名气更该主导选型。

方程求根：稳与快的搭档

换热器设计工单：污垢热阻随温度非线性变化，要解一个超越方程求稳态工作温度。二分法慢但绝对可靠（区间内有根就必然收敛）；牛顿法快（误差平方收敛——有效数字每步翻倍）但要求导数、且可能发散。生产代码的标准姿势是混合：先用二分稳住区间，每步尝试牛顿，牛顿结果跳出区间就退回二分。

import numpy as np

def design_eq(T, Ta=25.0, q0=8.0, k=0.05):
    """能量平衡方程：导热散热 = 内部发热（随温度增长）"""
    return k*(T - Ta) - q0*np.exp(-((T-140)/60)**2) + 0.0

def bisect(f, lo, hi, tol=1e-10):
    flo = f(lo)
    for _ in range(200):
        mid = 0.5*(lo + hi)
        if f(mid)*flo <= 0:
            hi = mid
        else:
            lo, flo = mid, f(mid)
        if hi - lo < tol:
            break
    return 0.5*(lo + hi)

def newton(f, x0, tol=1e-12, h=1e-6):
    x = x0
    for _ in range(60):
        fx = f(x)
        dfx = (f(x+h) - f(x-h)) / (2*h)     # 数值导数
        step = fx / dfx
        x_new = x - step
        if abs(x_new - x) < tol:
            return x_new
        x = x_new
    return x

root_b = bisect(design_eq, 30, 200)
root_n = newton(design_eq, 120.0)
print(f"二分解: {root_b:.8f}")
print(f"牛顿解: {root_n:.8f}")
print(f"方程残差: {abs(design_eq(root_n)):.2e}")
​

牛顿法的平方收敛值得眼见为实：把迭代过程打印出来，有效数字大致按 1、2、4、8 的节奏翻倍。但初值 120 换成 400 它可能一路飞出物理范围——快的代价是对起点的信任。混合策略（科学计算库求根器的内核做法）用二分区间当"安全绳"，牛顿当"加速器"，两者各取所长。

另一类翻车现场是对有重根或导数接近零的方程强推牛顿法：切线近乎水平，一步跳到天边。发现牛顿步长异常大时立刻回退到二分步，是健壮性代码的标配防御。

💡 关键直觉：误差阶是数值算法的"汇率"。同样是把区间加密一倍，二阶方法误差降四倍，四阶降十六倍——但高阶往往要求被积函数更光滑。光滑度换阶数，是这三类方法共同的交易结构。

三大常备药选型

本节要点回顾

• 龙格现象宣判等距高次插值死刑，多点场景交给分段样条；
• 带噪数据不插值，改用最小二乘拟合这条铁则呼应第 2 章；
• 误差阶是汇率：梯形二阶、辛普森四阶，光滑度是入场券；
• 自适应积分把工时自动集中到陡峭区，容差即预算；
• 整周期函数配梯形是隐藏的高精度通道，傅里叶计算的地基；
• 求根混合策略：二分当安全绳、牛顿当加速器，快而不断线。