教你如何利用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)同时也要留意容差不能设得比浮点精度还小。当容忍度小于机器精度时,迭代永远达不到条件,程序只会撞上最大迭代次数后失败。
常见坑点
- 直接用
==比较浮点数。0.1 + 0.2 == 0.3是False,任何基于相等判断的循环或分支都会失灵。
❌if x == 0.3:✅if math.isclose(x, 0.3, rel_tol=1e-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 的情况并报错,或换一个初值 / 改用二分法。
- 数值积分分段数给得太少就下结论。段数不足时截断误差很大,结果离真值可能差很远。
❌ 用n=10算完就当成精确值 ✅ 做"分段数翻倍"的收敛性检查,结果稳定后再采信。
- 把容差设得比机器精度还小。
tol=1e-20这种要求永远无法满足,循环必然撞上最大迭代次数。
❌tol=1e-18对双精度浮点做绝对误差判断 ✅ 容差与量级挂钩,用相对容差或取到1e-12量级即可。
- 忽略大数相减带来的抵消。两个接近的大数相减,有效数字几乎全丢,结果相对误差可能极大。
❌ 直接用sqrt(x**2 + 1) - x在 x 很大时计算 ✅ 改写成等价的稳定形式,例如先乘共轭表达式化简,避免相近量相减。
总结
| 任务 | 标准库方案 | 关键注意点 |
|---|
| 精确十进制 | decimal.Decimal | 必须用字符串构造 |
| 精确分数 | fractions.Fraction | 适合小分母的精确运算 |
| 复数 | cmath | 求负数的平方根要用它 |
| 求根 | 二分法 / 牛顿法 | 二分法要保证异号;牛顿法要防导数为零 |
| 积分 | 梯形法 / 辛普森法 | 用"分段翻倍"检查收敛性 |
| 比较 | math.isclose() | 不要用==比浮点 |
| 收敛判据 | 区间宽度或相对步长 | 别只看函数值 |
数值分析的核心不是记住公式,而是理解"近似有多准、误差从哪来"。写代码之前先想清楚这个问题的误差来源,比挑一个所谓"更高级"的算法重要得多;写完之后做一次收敛性验证,比多写几行优化更有价值。