前言数值分析numerical analysis研究的是用有限步的算术运算去近似求解数学问题求方程的根、算定积分的值、解线性方程组、对函数做插值。它和符号计算的分工很清楚——符号计算给的是精确表达式数值分析给的是一个足够接近的近似值和它的误差范围。新手最常见的误解有两个。第一是以为数值计算就等于浮点计算浮点计算就等于精确。事实上浮点数只是实数的离散近似0.1 0.2在二进制浮点下并不等于0.3而且每一步运算都会累积舍入误差。第二是以为标准库没有科学计算函数就做不了数值分析。math、decimal、fractions、cmath这几个标准库模块已经覆盖了大量基础数值任务本文的示例全部只用它们。本文按误差从哪里来怎么求根怎么求积分三条线展开给出可以直接运行并自己动手调参数的代码并说明每一处在数值上需要提防什么。一、误差从哪里来数值计算的误差主要有三个来源模型误差用简化模型近似真实现象、截断误差用有限项近似无穷过程比如把定积分换成有限个矩形面积之和、舍入误差浮点数位数有限每一步都要把结果舍入到可表示的最近值。判定一个近似值好不好不能只看差多少要看相对于量级差多少# 适用于 Python 3.8import matha 1e9b a 1e-3print(a b) # 可能为 True小数在大数上被吃掉print(math.isclose(a, b)) # 判断相对误差是否落在容忍范围内当两个数相差很大时把它们相减会发生抵消cancellation有效数字大量丢失相对误差被放大。数值分析里很多看起来正确的公式在浮点下失效根子就在抵消。误差类型来源缓解办法模型误差用简化模型代替真实系统靠领域知识评估程序帮不上忙截断误差有限项近似无穷过程增加项数 / 提高阶数直到结果稳定舍入误差浮点数位数有限避免大数相减、控制迭代次数抵消放大相近的大数相减改写公式或提高精度类型当需要精确的十进制结果时用decimal模块需要精确的分数运算时用fractions模块# 适用于 Python 3.8from decimal import Decimal, getcontextfrom fractions import Fractionprint(0.1 0.2) # 0.30000000000000004print(Decimal(0.1) Decimal(0.2)) # 0.3print(Fraction(1, 3) Fraction(1, 6)) # 1/2getcontext().prec 30 # 提高 Decimal 有效位数print(Decimal(1) / Decimal(3))注意Decimal的构造Decimal(0.1)用字符串才算精确写成Decimal(0.1)会先把二进制浮点转成十进制误差照样带进来。二、求方程的根求根最稳妥的方法是二分法bisection。前提是函数在区间两端取值异号中值定理保证区间内至少有一个根然后不断对半缩小区间。它收敛慢但几乎不会失败。# 适用于 Python 3.8def bisect(f, lo, hi, tol1e-10, max_iter200):flo f(lo)fhi f(hi)if flo 0:return loif fhi 0:return hiif flo * fhi 0:raise ValueError(区间两端必须异号或端点就是根)for _ in range(max_iter):mid (lo hi) / 2fmid f(mid)if fmid 0 or (hi - lo) / 2 tol:return midif flo * fmid 0:hi, fhi mid, fmidelse:lo, flo mid, fmidraise RuntimeError(达到最大迭代次数仍未收敛)print(bisect(lambda x: x * x - 2, 0, 2)) # 逼近 2 的平方根收敛判据用的是区间半宽(hi - lo) / 2小于容忍度tol而不是函数值——当函数在根附近非常平缓时函数值很小并不代表离根很近。比二分法快得多的是牛顿法Newtons method它用切线逼近迭代式是 x 减去 f(x) 除以 f(x)。但它需要导数且初值选得不好会发散或在导数接近 0 处爆炸。# 适用于 Python 3.8def newton(f, df, x0, tol1e-12, max_iter100):x x0for _ in range(max_iter):fx f(x)if abs(fx) tol:return xdfx df(x)if dfx 0:raise ZeroDivisionError(导数为 0牛顿法失效)x x - fx / dfxraise RuntimeError(未在迭代次数内收敛)print(newton(lambda x: x * x - 2, lambda x: 2 * x, 1.0))实践中的组合策略是先用二分法把区间缩到很小、拿到一个可靠初值再切换牛顿法快速收敛兼顾稳健和速度。三、数值积分定积分没有初等原函数时就要数值近似。最基础的是梯形法trapezoidal rule把积分区间切成 n 段每段用一个梯形面积近似再把它们加起来。# 适用于 Python 3.8def trapezoid(f, a, b, n1000):if n 0:raise ValueError(分段数必须为正)h (b - a) / ntotal (f(a) f(b)) / 2for i in range(1, n):total f(a i * h)return total * himport mathprint(trapezoid(math.sin, 0.0, math.pi, 1000)) # 逼近 2.0推演把 0 到 π 的 sin 曲线切成 1000 个等宽条端点函数值系数是 1/2中间各点系数是 1最后整体乘以宽度 h。因为 sin 在 [0, π] 上先增后减、形状接近对称梯形法在这里收敛得比较快换成在端点附近变化剧烈的函数同样的分段数误差会明显变大。更精确的是辛普森法Simpsons rule它用抛物线而不是直线去拟合每一小段相同分段数下误差通常更小——前提是分段数 n 为偶数。数值积分有一个必须养成的习惯把分段数翻倍看结果变化多少。如果翻倍后结果在小数点后第 6 位就稳定了说明当前精度够用如果还在跳说明分段不够或函数本身有奇异点。# 适用于 Python 3.8def adaptive_check(f, a, b):prev trapezoid(f, a, b, 100)for n in (200, 400, 800):cur trapezoid(f, a, b, n)print(fn{n:4} 结果{cur:.10f} 与上一次之差{abs(cur - prev):.2e})prev curadaptive_check(math.sin, 0.0, math.pi)同时也要留意容差不能设得比浮点精度还小。当容忍度小于机器精度时迭代永远达不到条件程序只会撞上最大迭代次数后失败。常见坑点直接用比较浮点数。0.1 0.2 0.3是False任何基于相等判断的循环或分支都会失灵。❌if x 0.3:✅if math.isclose(x, 0.3, rel_tol1e-9):。用Decimal(0.1)代替Decimal(0.1)。前者先把二进制浮点转成十进制误差已经被带进去了Decimal的精确性荡然无存。❌Decimal(0.1)并期待得到 0.1 精确值 ✅Decimal(0.1)用字符串构造。二分法不检查区间两端是否异号。两端同号时区间内不一定有根程序会一路收敛到一个假的结果且不报错。❌ 随便给个区间就调用求根函数 ✅ 先判断f(lo) * f(hi) 0否则给出明确错误。用函数值足够小当收敛判据。函数在根附近很平缓时离根还很远函数值就已经很小会提前误判收敛。❌if abs(f(x)) tol: return x作为唯一判据用于二分法时尤其危险 ✅ 二分法用区间宽度判据牛顿法可配合相对步长判据一起使用。牛顿法初值选在导数为 0 附近。除以接近零的导数会产生巨大跳步直接发散。❌ 对f(x) x**3取初值 0 后直接迭代 ✅ 检测导数为 0 的情况并报错或换一个初值 / 改用二分法。数值积分分段数给得太少就下结论。段数不足时截断误差很大结果离真值可能差很远。❌ 用n10算完就当成精确值 ✅ 做分段数翻倍的收敛性检查结果稳定后再采信。把容差设得比机器精度还小。tol1e-20这种要求永远无法满足循环必然撞上最大迭代次数。❌tol1e-18对双精度浮点做绝对误差判断 ✅ 容差与量级挂钩用相对容差或取到1e-12量级即可。忽略大数相减带来的抵消。两个接近的大数相减有效数字几乎全丢结果相对误差可能极大。❌ 直接用sqrt(x**2 1) - x在 x 很大时计算 ✅ 改写成等价的稳定形式例如先乘共轭表达式化简避免相近量相减。总结任务标准库方案关键注意点精确十进制decimal.Decimal必须用字符串构造精确分数fractions.Fraction适合小分母的精确运算复数cmath求负数的平方根要用它求根二分法 / 牛顿法二分法要保证异号牛顿法要防导数为零积分梯形法 / 辛普森法用分段翻倍检查收敛性比较math.isclose()不要用比浮点收敛判据区间宽度或相对步长别只看函数值数值分析的核心不是记住公式而是理解近似有多准、误差从哪来。写代码之前先想清楚这个问题的误差来源比挑一个所谓更高级的算法重要得多写完之后做一次收敛性验证比多写几行优化更有价值。