教你如何利用python进行数值分析

前言

数值分析（numerical analysis）研究的是"用有限步的算术运算去近似求解数学问题"：求方程的根、算定积分的值、解线性方程组、对函数做插值。它和符号计算的分工很清楚——符号计算给的是精确表达式，数值分析给的是一个足够接近的近似值和它的误差范围。

新手最常见的误解有两个。第一是以为"数值计算就等于浮点计算，浮点计算就等于精确"。事实上浮点数只是实数的离散近似，0.1 + 0.2在二进制浮点下并不等于0.3，而且每一步运算都会累积舍入误差。第二是以为"标准库没有科学计算函数就做不了数值分析"。math、decimal、fractions、cmath这几个标准库模块已经覆盖了大量基础数值任务，本文的示例全部只用它们。

本文按"误差从哪里来""怎么求根""怎么求积分"三条线展开，给出可以直接运行并自己动手调参数的代码，并说明每一处在数值上需要提防什么。

一、误差从哪里来

数值计算的误差主要有三个来源：模型误差（用简化模型近似真实现象）、截断误差（用有限项近似无穷过程，比如把定积分换成有限个矩形面积之和）、舍入误差（浮点数位数有限，每一步都要把结果舍入到可表示的最近值）。

判定一个近似值好不好，不能只看"差多少"，要看相对于量级差多少：

# 适用于 Python 3.8+
import math

a = 1e9
b = a + 1e-3
print(a == b) # 可能为 True：小数在大数上被"吃掉"
print(math.isclose(a, b)) # 判断相对误差是否落在容忍范围内

当两个数相差很大时，把它们相减会发生抵消（cancellation）：有效数字大量丢失，相对误差被放大。数值分析里很多"看起来正确的公式"在浮点下失效，根子就在抵消。

当需要精确的十进制结果时，用decimal模块；需要精确的分数运算时，用fractions模块：

# 适用于 Python 3.8+
from decimal import Decimal, getcontext
from fractions import Fraction

print(0.1 + 0.2) # 0.30000000000000004
print(Decimal("0.1") + Decimal("0.2")) # 0.3
print(Fraction(1, 3) + Fraction(1, 6)) # 1/2

getcontext().prec = 30 # 提高 Decimal 有效位数
print(Decimal(1) / Decimal(3))

注意Decimal的构造：Decimal("0.1")用字符串才算精确，写成Decimal(0.1)会先把二进制浮点转成十进制，误差照样带进来。

二、求方程的根

求根最稳妥的方法是二分法（bisection）。前提是函数在区间两端取值异号（中值定理保证区间内至少有一个根），然后不断对半缩小区间。它收敛慢但几乎不会失败。

# 适用于 Python 3.8+
def bisect(f, lo, hi, tol=1e-10, max_iter=200):
flo = f(lo)
fhi = f(hi)
if flo == 0:
return lo
if fhi == 0:
return hi
if flo * fhi > 0:
raise ValueError("区间两端必须异号，或端点就是根")
for _ in range(max_iter):
mid = (lo + hi) / 2
fmid = f(mid)
if fmid == 0 or (hi - lo) / 2 < tol:
return mid
if flo * fmid < 0:
hi, fhi = mid, fmid
else:
lo, flo = mid, fmid
raise RuntimeError("达到最大迭代次数仍未收敛")

print(bisect(lambda x: x * x - 2, 0, 2)) # 逼近 2 的平方根

收敛判据用的是区间半宽(hi - lo) / 2小于容忍度tol，而不是函数值——当函数在根附近非常平缓时，"函数值很小"并不代表"离根很近"。

比二分法快得多的是牛顿法（Newton's method），它用切线逼近，迭代式是 x 减去 f(x) 除以 f'(x)。但它需要导数，且初值选得不好会发散或在导数接近 0 处爆炸。

# 适用于 Python 3.8+
def newton(f, df, x0, tol=1e-12, max_iter=100):
x = x0
for _ in range(max_iter):
fx = f(x)
if abs(fx) < tol:
return x
dfx = df(x)
if dfx == 0:
raise ZeroDivisionError("导数为 0，牛顿法失效")
x = x - fx / dfx
raise RuntimeError("未在迭代次数内收敛")

print(newton(lambda x: x * x - 2, lambda x: 2 * x, 1.0))

实践中的组合策略是：先用二分法把区间缩到很小、拿到一个可靠初值，再切换牛顿法快速收敛，兼顾稳健和速度。

三、数值积分

定积分没有初等原函数时，就要数值近似。最基础的是梯形法（trapezoidal rule）：把积分区间切成 n 段，每段用一个梯形面积近似，再把它们加起来。

# 适用于 Python 3.8+
def trapezoid(f, a, b, n=1000):
if n <= 0:
raise ValueError("分段数必须为正")
h = (b - a) / n
total = (f(a) + f(b)) / 2
for i in range(1, n):
total += f(a + i * h)
return total * h

import math
print(trapezoid(math.sin, 0.0, math.pi, 1000)) # 逼近 2.0

推演：把 0 到 π 的 sin 曲线切成 1000 个等宽条，端点函数值系数是 1/2，中间各点系数是 1，最后整体乘以宽度 h。因为 sin 在 [0, π] 上先增后减、形状接近对称，梯形法在这里收敛得比较快；换成在端点附近变化剧烈的函数，同样的分段数误差会明显变大。

更精确的是辛普森法（Simpson's rule），它用抛物线而不是直线去拟合每一小段，相同分段数下误差通常更小——前提是分段数 n 为偶数。

数值积分有一个必须养成的习惯：把分段数翻倍，看结果变化多少。如果翻倍后结果在小数点后第 6 位就稳定了，说明当前精度够用；如果还在跳，说明分段不够或函数本身有奇异点。

# 适用于 Python 3.8+
def adaptive_check(f, a, b):
prev = trapezoid(f, a, b, 100)
for n in (200, 400, 800):
cur = trapezoid(f, a, b, n)
print(f"n={n:>4} 结果={cur:.10f} 与上一次之差={abs(cur - prev):.2e}")
prev = cur

adaptive_check(math.sin, 0.0, math.pi)

同时也要留意容差不能设得比浮点精度还小。当容忍度小于机器精度时，迭代永远达不到条件，程序只会撞上最大迭代次数后失败。

常见坑点

1. 
   
2. 直接用==比较浮点数。0.1 + 0.2 == 0.3是False，任何基于相等判断的循环或分支都会失灵。
3. 
   

❌if x == 0.3:✅if math.isclose(x, 0.3, rel_tol=1e-9):。

1. 
   
2. 用Decimal(0.1)代替Decimal("0.1")。前者先把二进制浮点转成十进制，误差已经被带进去了，Decimal的精确性荡然无存。
3. 
   

❌Decimal(0.1)并期待得到 0.1 精确值 ✅Decimal("0.1")，用字符串构造。

1. 
   
2. 二分法不检查区间两端是否异号。两端同号时区间内不一定有根，程序会一路收敛到一个假的结果，且不报错。
3. 
   

❌ 随便给个区间就调用求根函数 ✅ 先判断f(lo) * f(hi) <= 0，否则给出明确错误。

1. 
   
2. 用"函数值足够小"当收敛判据。函数在根附近很平缓时，离根还很远函数值就已经很小，会提前误判收敛。
3. 
   

❌if abs(f(x)) < tol: return x作为唯一判据（用于二分法时尤其危险） ✅ 二分法用区间宽度判据，牛顿法可配合相对步长判据一起使用。

1. 
   
2. 牛顿法初值选在导数为 0 附近。除以接近零的导数会产生巨大跳步，直接发散。
3. 
   

❌ 对f(x) = x**3取初值 0 后直接迭代 ✅ 检测导数为 0 的情况并报错，或换一个初值 / 改用二分法。

1. 
   
2. 数值积分分段数给得太少就下结论。段数不足时截断误差很大，结果离真值可能差很远。
3. 
   

❌ 用n=10算完就当成精确值 ✅ 做"分段数翻倍"的收敛性检查，结果稳定后再采信。

1. 
   
2. 把容差设得比机器精度还小。tol=1e-20这种要求永远无法满足，循环必然撞上最大迭代次数。
3. 
   

❌tol=1e-18对双精度浮点做绝对误差判断 ✅ 容差与量级挂钩，用相对容差或取到1e-12量级即可。

1. 
   
2. 忽略大数相减带来的抵消。两个接近的大数相减，有效数字几乎全丢，结果相对误差可能极大。
3. 
   

❌ 直接用sqrt(x**2 + 1) - x在 x 很大时计算 ✅ 改写成等价的稳定形式，例如先乘共轭表达式化简，避免相近量相减。

总结

数值分析的核心不是记住公式，而是理解"近似有多准、误差从哪来"。写代码之前先想清楚这个问题的误差来源，比挑一个所谓"更高级"的算法重要得多；写完之后做一次收敛性验证，比多写几行优化更有价值。