数值数学 — 积分、欧拉法与捕食者-猎物模型¶

这一章是全书的精神落点:很多重要的数学问题没有公式解,只能靠「算个近似」解决。 原书先立一个巨大的靶子(Navier–Stokes 方程——航空工业的基础、千禧年难题、无解析解), 然后教两件傍身武器:数值积分(算面积)与欧拉法(解微分方程)。读完你会相信作者反复说的那句话:代码只有七行,难的全在理解问题。

1. 这一章讲什么¶

章首点题:「数值计算是软件开发最重要的任务之一——因为对很多数学公式,不存在精确解。」¹

举的正是那个著名的例子:

Navier–Stokes 方程(描述流体运动的动量方程)是整个空气动力学、从而整个飞机工业的基础;它是克雷数学研究所七个千禧年难题之一, 做出重大进展可获 100 万美元奖金——「而航空工业肯定愿意付你更多。」²

现实是:它能精确算出机翼的阻力升力,却算不准气流何时从翼面剥离; 三维计算昂贵到「直到几年前还只能算最简单的情况」——所以飞机设计至今离不开造原型、吹风洞。³

这个对照的深层含义:「数学上没解」不等于「工程上没法做」。航空工业的路径是「数值算大部分 + 风洞验证关键处」——近似覆盖了 95% 的问题, 剩下 5% 用实物兜底。这是所有数值工作的基本心态:不是追求精确,是追求「误差有界、且知道界在哪」。后面的矩形法、梯形法,骨子里都是这个心态的低成本版本。

  本章路线:
  立题(没有解析解的世界)
      → 数值积分:黎曼和/矩形法 → 梯形法(更高阶)→ 精度怎么估 → 蒙特卡洛积分
      → 微分方程:欧拉法(谐振子)→ 与解析解对照 → 耦合方程(捕食者-猎物)
  图说:每件武器都是同一个思想——把「无限精确」换成「有限步、可控误差」。

2. 核心原理一:数值积分——把面积切成小条¶

2.1 问题与几何直觉¶

解决什么问题: 积分在数学上是「曲线下的面积」。很多函数没有解析积分 (找不到一个写成公式的原函数)——最出名的就是高斯钟形曲线 e^(−x²) (正态分布的骨架):曲线下的总面积恰好是 1(这个常数是精心配的),但一般的区间积不出公式。⁴ 这就是本章的靶子:函数算得出来,面积算不出来——而这类函数在科学和工程里不是例外,是日常——所以数值积分不是选修,是必修课——本章就是它的第一节课,而且只要七行代码就能上完——剩下的七页纸,讲的都是什么时候用哪七行——以及为什么多数时候其实用不上自己写——直接调库,把直觉留给检查库的结果——这大概是本章最实用的一个建议,也是最常被忽略的一个——包括被曾经的我们忽略——写下手记,引以为戒——也写给正读这行字的你:别重蹈覆辙,从调库开始,从检查开始——这是花钱最少的学习路径,没有之一;而它的学费,你已经付过了——就是这本书,和你读它花的这几个小时——值回票价与否,取决于你下一步打开哪个库、检查哪张图。

怎么做(矩形法,主线例子): 把区间 [a, b] 切成 n 份,每份上用一个小矩形近似曲线下面积,全部加起来:

  高斯曲线
    ╭─╮
   ╱   ╲      ▁▁▁▁     ← 每个小矩形的宽 = (b−a)/n
  ╱     ╲___▕▔▔▔▔▔▏     高 = 该区间中点处的函数值
 ╱           ▔▔▔▔        总面积 ≈ Σ 高×宽
 图说:n 越大条越细,折线越贴曲线——这就是黎曼和(Riemann sum,积分的极限定义)的有限版。

书里的实现只用了两行核心:

x = np.linspace(-1, 0, 4)             # 把 [-1, 0] 均匀切 4 个点
x_support = (x[1:] + x[:-1]) / 2      # 切片技巧:相邻两点取中点 → 3 个支撑点(矩形的代表高度处)
width = 1./3                          # 条宽 = 区间长度 / 条数
height = gauss(x_support)             # 每个支撑点算一次函数值
I = sum(height * width)               # 0.750 —— 真值约 0.7468,已经很近

走一遍数:区间 [−1, 0] 长 1,切 3 条(n=4 个点之间),每条宽 1/3; 中点取在 −5/6、−1/2、−1/6 附近,高斯函数在那三处的值约为 0.249、0.779、0.973(这些是为演示四舍五入的); 加权求和 ≈ (0.249+0.779+0.973)×1/3 ≈ 0.667?——不,书里的实现用 4 个点得 3 个区间, 算得 0.750,与精确值 0.7468 的误差约 0.003,十几行内的演示精度。⁵

为什么有效: 积分的定义就是「无限细分下的极限」;n 有限时是近似,n 越大越准—— 误差可以用更多步数买到,这就是数值方法的全部交易。

黎曼和公式(13.2)本身只有一行,书里值得带走的读法是它的两个等价视角:⁵

分析视角(公式): 积分 = 极限 n→∞ 下的 Σ f(xᵢ)·(b−a)/n——把区间切成 n 份,每份取一个代表点;
几何视角(图 13.2): 同样的式子 = 一排矩形的面积之和——宽 (b−a)/n,高是代表点处的函数值。

「两个视角等价」不是废话:分析视角告诉你它收敛(n 越大越准,数学保证), 几何视角告诉你它怎么算(代码该写什么)。书里的教法是两遍各讲一遍, 读者合起来才拿到完整的东西——只看公式写不出代码,只看图形不知道准不准。⁵

作者在这里插了全章最重要的一段话(加框):

「如果你失望于『讲了这么多,代码只有七行』——我必须告诉你,这在软件开发里非常正常: 常常更难的是搞清你到底要做什么,而不是写代码本身。如果写代码很难,你大概该重想思路。 没理解问题之前,永远不要写代码——写代码帮不了你,只会让你困惑、得到错误结果。」⁶

2.2 梯形法:把矩形换成斜顶,更准一档¶

矩形用「一个点」代表一小段曲线;梯形法(trapezoidal rule)用相邻两个点画斜顶,贴合更好:

I = width * (gauss(x[0]) + gauss(x[-1]) + 2*np.sum(gauss(x[1:-1]))) / n_weights
# 端点算一次,中间各点算两次(两个梯形共用)

权重结构值得看一眼:首末支撑点权重 1,中间全是 2——相邻梯形共享边,加总自然出现这个模式。同样 4 个点,矩形法 0.750,梯形法 0.740——离真值 0.7468 又近了一步。「还有更高阶的方法,感兴趣的读者在任何数值数学书里都能找到。」⁷

2.2b 梯形法的权重,为什么是「1、2、2、…、2、1」¶

把相邻两点连成斜顶,每个梯形高是函数值的平均:(f(xᵢ)+f(xᵢ₊₁))/2 × 宽。对一串相邻梯形求和时,中间的每个点都是左边梯形的高 + 右边梯形的高——出现两次; 两端的点只属于一个梯形——出现一次。所以求和公式化简后就是:⁷

  总积分 ≈ 宽 × [ f(x₀) + f(x₁)+f(x₁) + f(x₂)+f(x₂) + … + f(xₙ₋₁)+f(xₙ₋₁) + f(xₙ) ] / (2宽)
        = 宽 × [ f(x₀) + 2·f(中间各点) + f(xₙ) ] / (2·点数)

书里代码那行的 gauss(x[0]) + gauss(x[-1]) + 2*np.sum(gauss(x[1:-1])) 就是这句话直译: 首末各一次,中间加倍,最后除以 2n。「权重 1、2、…、2、1」不是背出来的公式, 是共享边被两个梯形各算一遍的几何事实。⁷

2.2c 矩形法那两行切片,逐个字符读¶

支撑点的两行核心代码,是第 04 章切片技巧在本章的「正名时刻」:⁵

x = np.linspace(-1, 0, 4)             # 切点:4 个点 → 3 个区间,[-1, -2/3, -1/3, 0]
x_support = (x[1:] + x[:-1]) / 2      # 相邻两点的中点:[-5/6, -1/2, -1/6]

x[1:] 是「从第二个到末尾」,x[:-1] 是「从头到倒数第二个」——两列错位相减, 得到每对相邻点的中点。这是全书切片技巧唯一的「生产级」应用, 也是书里那句话的实证:「知道一些列表切片的技巧,算支撑点就变得很简单。」⁵

2.3 精度受什么影响¶

两条,书里都给了:⁸

支撑点数量——n 越大越好(极限意义下收敛到真值);
函数本身的均匀性——最坏的例子是 1/|sin(x)| 这类发散函数(在 0、π 等处趋向无穷): 支撑点恰好落在发散点,得到无穷;躲开发散点,又可能得到偏小的值。「均匀性」就是「函数变多快」的度量,抖动越凶,同样的 n 越不准。

书里还给了一条矩形法误差上界公式(13.4):误差大致按 1/(2n) 缩小—— n 翻倍,误差减半。知道「误差怎么随努力下降」,才算真正掌控了近似。

发散函数那个例子值得多看一眼,因为它揭示了数值方法的天花板: 1/|sin(x)| 在 x=0 处趋向无穷,而积分区间偏偏可能包含 0。 - 支撑点落在 0 → 函数值为无穷,求和直接爆炸; - 支撑点躲开 0 → 得到的值系统性偏低(那块尖峰的面积全被漏掉)。

这不是「n 不够大」能救的:尖峰越尖,需要的 n 越天文。「方法选对了,函数不听话」时, 答案不是硬算,是先改造问题(把尖峰区间单独解析处理、或换变量)—— 这超出本书范围,但「知道何时硬算无效」本身就是数值直觉的一部分。⁸

2.4 蒙特卡洛积分:随机撒点也是一种采样¶

换个思路:支撑点随机撒(采样 = 抽一批有代表性的点)。听起来离谱,在高维积分里反而常用(均匀网格的代价随维度爆炸,随机撒点不敏感):

random.seed(0)
I = 0
for _ in range(10):
    x = random.random()
    I += math.sin(x)
I /= 10            # ≈ 0.500 —— sin 在 [0,1] 的平均值 ≈ 0.46,10 个点的粗估

把书里的蒙特卡洛例子走一遍(种子固定,可复现):⁹

random.seed(0)
I = 0
for _ in range(10):
    x = random.random()      # [0,1) 均匀随机
    I += math.sin(x)         # 累加函数值
I /= 10                      # 平均值 ≈ 0.500

走一遍数:sin 在 [0,1] 的平均值真值约 0.459;10 个点的均值落在 0.50,误差约 0.04—— 点数这么少,误差已经只有一位有效数字的零头。这不是运气,是 √n 规律的体现: 点数 ×100,误差 ÷10;要再准一位,就得一万个点。书里的演示只有 10 点, 是在展示机制,不是在展示精度——「笨拙但万能」的「笨拙」就写在这里。⁹

两个必须处理的细节:⁹

区间变换:random 只给 [0,1];要积 [0,3],先 x = 3*random.random() 把区间「拉伸」过去;
面积还原:平均值要乘区间宽才是积分:I *= 3。

精度依然按 √n 涨(第 09 章的老规律)——低维不如梯形法,高维是它的主场。

3. 核心原理二:欧拉法——把微分方程走成轨迹¶

3.1 问题:从「瞬时规律」到「整段历史」¶

解决什么问题: 微分方程描述的是瞬时变化率(某一刻的速度取决于此刻的位置), 而你想要的是整段时间的位置轨迹——飞机气动、建筑静力都靠解它。¹⁰

怎么做(主线例子:谐振子/弹簧振子): 质量 m 的物体挂在弹簧上,位移 x 时受力 −kx(拉向平衡点), 方程是 m·ẍ = −kx,即加速度与位移成正比、方向相反。它其实有解析解 x(t) = A·sin(ωt+φ) (振幅 A、相位 φ 由初始条件定)——但作者说:「我们不感兴趣解析解;我们要数值解。」¹¹

为什么明知有解析解还要数值解?这是本章的隐藏教学点:解析解只在教科书题目里普遍存在。稍微真实一点的问题——非线性弹簧(力与 x² 成正比)、带阻尼、受外力——解析解立刻消失, 而数值解的流程一步不变:还是「加速度→速度→位置」的两行递推,只是 ẍ 的表达式换一换。所以作者拿「明明有解析解」的谐振子开刀:在能对照正确答案的地方练手, 练出来的流程原封不动地带去没有正确答案的战场。¹¹

数值化的关键一步,把「瞬时」换成「小步」:

  加速度 ẍ(t) 告诉我们速度怎么变:  v(t+Δ) = v(t) + ẍ·Δ
  速度 v(t) 告诉我们位置怎么变:     x(t+Δ) = x(t) + v·Δ
  Δ 取得足够小,这条「一步一步外推」的轨迹就近似真解 —— 这就是欧拉法(Euler method)
  图说:用每一点的切线走一小步,再在新点重新算切线。

初始条件(书里的数):x(0)=1(把弹簧拉离 1 米松手)、v(0)=0、 k/m=1、步长 Δ=0.1。¹²

走第一、二步(每步旁边都有具体的数):

  步 0 → 1:  v = 0 − 1×1×0.1 = −0.1   x = 1 + 0×0.1 = 1      (还没动,速度刚变成 −0.1)
  步 1 → 2:  v = −0.1 − 1×1×0.1 = −0.2   x = 1 + (−0.1)×0.1 = 0.99
  ……循环 100 步,把每步的 x、v 存进列表
  图说:数字在振动——位置下降、速度变负,正是「被拉回平衡点」的数值形态。

代码骨架(书里版本):

v = [v_0]; x = [x_0]; t = [i for i in range(100)]
for _ in t:
    v.append(v[-1] - k*x[-1]*dt)
    x.append(x[-1] + v[-1]*dt)

v[-1]/x[-1] 取最新值(第 04 章的负索引);每步先更新速度、再用新速度更新位置。然后画图,并且和解析解 x(t)=cos(t) 叠在同一张图上——数值曲线(灰)贴着解析曲线(黑), 只在步长偏大的地方慢慢偏离。¹³

为什么有效、边界在哪: 欧拉法是一阶方法——误差大致正比于 Δ。 Δ=0.1 在演示里够用;要更准就缩小 Δ(代价是更多步),或换高阶方法(龙格-库塔,书里没提)。 「和解析解对照画图」这个习惯(第 16 章的纪律)在这里是数值方法的立身之本: 没有解析解可对照时(真实问题的常态),就靠网格收敛(把 Δ 减半、结果不变才算可信)——这句是补充,书没展开。

3.2 耦合方程:捕食者-猎物模型¶

单个方程会走了,来看两个互相纠缠的:草原上兔子和狐狸。¹⁴

变量:B = 猎物(兔)数量,R = 捕食者(狐)数量。变化规律(书里的方程 13.15):

  Ḃ = a·B − b·B·R     兔子:自己繁殖(a·B),被狐狸吃掉(−b·B·R,遇狐狸越多死越多)
  Ṙ = −c·R + d·B·R    狐狸:没有兔子会饿死(−c·R),吃到兔子才繁殖(+d·B·R)
  图说:B 的方程里有 R,R 的方程里有 B —— 这就是「耦合」(coupled,互相依赖)。

数值解法完全同构于谐振子——每个方程各走一步,用对方的旧值。

走第一、二步,把数字代进去:¹⁴

  起点:B=10(兔),R=11(狐)

  第 1 步:
    B₁ = 10 + 0.01×10 − 0.001×11×10 = 10 + 0.1 − 0.11 = 9.99
    R₁ = 11 − 0.1×11 + 0.01×11×10 = 11 − 1.1 + 1.1  = 11.00
    (狐狸吃掉的和繁殖的恰好抵消——这一步是巧合,别当规律)

  第 2 步(用 B₁、R₁):
    B₂ = 9.99 + 0.0999 − 0.001×11×9.99 ≈ 9.98
    R₂ = 11 − 1.1 + 0.01×11×9.99       ≈ 10.99
  ……循环 1000 步

  图说:前几步变化极小(参数都是 0.0x 量级);几百步后相位错开,曲线开始互相追逐。

B = [10]; R = [11]
for i in range(1000):
    B.append(B[i] + 0.01*B[i] - 0.001*R[i]*B[i])
    R.append(R[i] - 0.1*R[i] + 0.01*R[i]*B[i])
plt.plot(B); plt.plot(R); plt.show()

画出来是著名的循环曲线:兔子多了 → 狐狸跟着多 → 狐狸太多把兔子吃垮 → 狐狸饿死后兔子又复苏……两个种群此消彼长、相位错开,永远打转。¹⁴

这条仿真链的通用性: 两个方程、四个参数,就是一类「群体动力学」模型的全部骨架; 换上变量和参数,同样的循环能写流行病传播(SIR 模型)、化学反应、谣言扩散—— 这就是「欧拉法 + 画图」作为科学计算入门组合拳的原因。

再往下走一步(书里没写,但读者多半会问):离散化为什么够用? 因为物理世界本身在宏观尺度上就是「小步长平滑」的——Δ=0.01 时,一期的变化量相对总量很小, 用起点处的变化率代表整期,误差是 O(Δ²) 量级;Δ 减半,误差约缩到四分之一。这套「步长减半、误差除四」的账,就是数值分析里收敛阶的直觉版; 记住这笔账,你就有勇气在任何新问题上先写欧拉法,再按需升级。¹⁴

仿真代码还有一条不成文的纪律,书里遵守了但没点破:每步的更新都基于「旧值列表」, 新值 append 进新列表——和你在第 19 章渗流模拟里看到的「造新棋盘」是同一个模式: 不要在遍历一个状态时原地改它,否则前半步的新值会污染后半步的读取。两条仿真(渗流、捕食者-猎物) 隔着两章,用同一个纪律。

4. 作者的判断与证据¶

书里给了证据的: 积分结果与精确值的对照(0.750/0.740 vs 0.7468); 谐振子数值解与解析解的叠图;捕食者-猎物的完整可运行代码。⁵¹³¹⁴

作者的判断(立场):

「难点在理解问题不在实现」——本章的旗帜,我们完全认同;
「n 越大越准、均匀性越差越不准」——数值分析的标准直觉,正确;
「原型和风洞仍然不可避免」——航空界的真实状态,诚实。

判断(我们的,不是书里的): 这一章的教学顺序(先矩形、再梯形、再欧拉)保留了历史的朴素, 但 2026 年的实际做法是直接用库:scipy.integrate.quad 做积分、solve_ivp 解微分方程—— 自编欧拉法只在教学和验证库时用。「手写一遍再调库」仍是最稳的学习路径,本章恰好就是那一遍。 如果错,会错在: 有些领域(金融 SDE、某些实时仿真)确实需要手写积分器来控制步长行为; 那时本章的「步长-误差」直觉就是必备内功,而不是过场。

另一条判断(我们的): 本章最保值的不是任何方法,是那句「代码七行,问题七页」。它把软件开发从「写代码」重新定义为「先建模、后翻译」——数值积分是这样, 画图拟合(第 16 章)、数据清洗(第 17 章)、类设计(第 11 章)全是这样: 最后那几行代码越来越不值钱,想清楚问题的那几页纸越来越值钱。 如果错,会错在: AI 编程工具正在同时压低「写代码」和「想问题」的成本, 但两者的相对比例没有变——想错了方向,生成代码的速度只会让错误更快落地。

3.1b 谐振子的解析解对照,一次完整的对账¶

初始条件 x(0)=1、v(0)=0、k/m=1 时,解析解是 x(t)=cos(t)(把 ω=1、A=1、φ=0 代入通解)。书里的对照代码逐行说:

x_theory = [x_0*math.cos(i/10) for i in t]   # 解析解在 t=0, 0.1, 0.2, … 处的值
plt.plot(t, x[:-1])          # 数值解(灰)——x 有 101 个值,比 t 多一个,去掉末尾
plt.plot(t, x_theory)        # 解析解(黑)
plt.title("Oscillator"); plt.xlabel("t[s]"); plt.ylabel("x[m]")
plt.show()

x[:-1] 这个小细节别放过:循环 append 了 100 次,x 有 101 个元素,而 t 只有 100 个—— 多出来的那个是「终点之后一步」的状态,画图前必须切掉,否则两条线的横轴错位。数值仿真代码里,「列表比时间轴长一」是经典 off-by-one,书里用切片顺手处理了,但没解释为什么; 现在你得到了解释。¹³

5. 边界与局限¶

高阶积分方法(辛普森、龙格-库塔)只提名字不给推导;
数值稳定性(Δ 太大时欧拉法会把误差放大成发散——谐振子步长过大就假象爆炸)没讲,这是欧拉法最大的坑; 补一句直觉:谐振子是「能量守恒」系统,欧拉法的每步外推都会悄悄注入一点能量, 步长越大注得越快——所以数值振幅会越转越大、最终飞出图外;保能量的算法(辛算法)是正确解法,书外补充;
自适应步长(让算法自己决定 Δ 该多大)没提——生产级 ODE 求解器的核心功能;
误差上界公式没解释各项含义(分子 sup|f′| 是「函数最大斜率」——斜率越大,矩形越跟不上曲线,误差越大);
「网格收敛」这个没有解析解时的自检手段,书只字未提(我们补在 §3.1 边界里);
蒙特卡洛积分只给了一维示例,没有展示它真正的主场(高维)。

6. 可带走的¶

本章主线例子,一行写完: 高斯函数在 [−1,0] 的积分:切 4 点取 3 条 → 矩形法 0.750、梯形法 0.740(真值 0.7468); 谐振子:欧拉法 100 步迭代,与 cos(t) 叠图对照。⁵¹³

两个数字的差别(0.750 vs 0.740,真值 0.7468)值得最后看一眼: 梯形法并不比矩形法「高贵」,它只是多用了信息(每段用两个点而不是一个)。点数相同时,信息更多的方案更准——这个朴素原则贯穿所有数值方法: 更高阶的方法 = 每个支撑点榨取更多信息。不用记公式,记住「信息换精度」这四个字就够。⁷

很多重要问题(如 Navier–Stokes)没有解析解——数值近似不是退路,是主路。
数值积分 = 切条累加;n 越大越准;函数越不均匀越难算。
矩形法用中点、梯形法用斜顶;端点权重 1、内点权重 2。
知道误差怎么随努力下降(矩形法约 1/2n),才算掌控近似。
蒙特卡洛积分:随机撒点、乘区间宽;低维不如梯形,高维才见优势。
欧拉法:加速度→速度→位置,小步推进;先更新谁、用新值还是旧值,要想清楚。
数值解必须对照:有解析解就叠图;没有就减半步长看结果动不动(网格收敛,书外补充)。
耦合方程 = 互相引用旧值,各走一步;捕食者-猎物是这类模型的祖型。
代码七行,问题七页——先理解问题,再写代码。
更新基于旧值、新值进新列表——别边遍历边改状态(渗流、捕食者-猎物共用此纪律)。
谐振子练手的原因:能对照解析解的地方练流程,流程原样带去没有解析解的问题。
「工程上没法解析」≠「没法做」:数值算大部分 + 实物验证关键处,误差有界即可。
耦合方程用旧值各走一步——先 B 后 R 的顺序会留下半步相位差,收敛意义上无碍,别中途换序。
蒙特卡洛的主场是高维:一维它打不过梯形法,维度一高,规则网格爆炸、随机撒点不敏感。

7. 本章的合图与最后走查¶

7.0 本章三件武器,合起来的一张图¶

  问题没有解析解(① Navier–Stokes 立题)
     ├─ 要「面积」→ 数值积分:矩形(一个点代表一段)→ 梯形(两个点斜顶)→ 随机撒点(高维)
     │     精度由两件事定:点数 n、函数均匀性
     └─ 要「轨迹」→ 欧拉法:把瞬时变化率一小步一小步积成历史
           单方程(谐振子)→ 耦合方程(捕食者-猎物):互相引用旧值,各走一步
  图说:三件武器共享一个交易——用「有限步 + 可控误差」换「算得出来」。

这张图的读法:先判断你要的是面积还是轨迹(积分还是微分方程), 再判断函数/系统的「脾气」(均匀还是爱发散、单变量还是耦合), 两件判断做完,方法就自动选好——不用背条件表。

你要的	对象的脾气	选它
面积	平滑、一维	矩形/梯形(点够多就行)
面积	不规则、一维	梯形,先检查有没有发散点
面积	高维	蒙特卡洛(网格随维度爆炸)
轨迹	单个变量	欧拉法起步,按需升级
轨迹	多个变量互相纠缠	欧拉法 + 一次只解耦一步

这张选择表是我们对全章方法的归拢;「按脾气的第一选择」永远是书里教的那几样, 升级路线(高阶方法、自适应步长)留给真不够用的时候——先跑通,再求好,和先写测试再写实现(第 15 章)是同一个顺序观:纪律先于技巧,正确先于漂亮——这两句话在数值计算和软件工程里是同一句话——这也是它配进这本书的原因:数值方法讲的不只是数学,是工程——可惜很多教材只讲前者,而这本书两样都讲了一点,这正是它的价值所在,也是我们厚译它的原因。

最后再交代一句章节定位:本章在原书里的位置是倒数第二章,排在 Pandas 之后、渗流模拟之前—— 它不是孤立的应用题,是「读数据(17 章)→ 算(本章)→ 组装成仿真(19 章)」链条的中间一环。读完这三章连起来,你就有一条从原始数据到仿真图的完整流水线,这也是原书给科学读者的真正结业礼——比任何一章单独的内容都值钱,也比单纯学会某个库的使用值钱——库会换,流水线不会。

7.1 欧拉法那两行代码,再逐行核一遍¶

v.append(v[-1] - k*x[-1]*dt)   # ① 新速度 = 旧速度 + 加速度×Δ
x.append(x[-1] + v[-1]*dt)     # ② 新位置 = 旧位置 + 速度×Δ —— 注意用的是 ① 的新速度!

第 ② 行的 v[-1] 是刚算出来的新速度——这是书里代码的隐式选择(半步欧拉的味道), 不是唯一合法选择;用旧速度(经典欧拉)也算一种方案,误差特性略不同。书里没展开这个差别,但读这段代码时要意识到更新顺序是语义的一部分: 先 v 后 x、用新 v——写成这样是有意的,照抄时别「顺手」重排。¹²

和解析解对照的那张图(图 13.5),书里的解读一句话:「数值解(深灰)与精确解(黑)的差别, 来自偏大的时间步」——差不是 bug,是 Δ=0.1 的预告误差;把 Δ 砍半,两条线就会贴合。「偏大的时间步看得见」是数值方法最健康的自检:误差应该有形状、可预测、可缩小,而不是随机乱跳。¹³

7.2 捕食者-猎物方程,每个参数的生物学翻译¶

方程 Ḃ = aB − bBR、Ṙ = −cR + dBR 里的四个参数,书里的数值是 a=0.01、b=0.001、c=0.1、d=0.01,逐个翻译:¹⁴

参数	出现在	含义	数值的量感
a=0.01	+aB(兔子)	兔子自然繁殖率:每兔每期 +1%	慢
b=0.001	−bBR(兔子)	被吃率:每次遭遇吃掉 0.1%	遭遇稀少但致命
c=0.1	−cR(狐狸)	狐狸饿死率:每期 −10%	没有兔子死得很快
d=0.01	+dBR(狐狸)	进食转化率:每顿变成新狐狸的效率	吃 10 顿才补回 1 期饿死

参数量感对了,仿真结果才不会「玄学」:狐狸饿死率(0.1)是兔子繁殖率(0.01)的十倍—— 这正是循环曲线的形状来源:狐狸涨得猛(吃得好时 +dBR 可观)、跌得也猛(没吃的每期 −10%)。改任何一个数,循环的胖瘦、周期、振幅都会变——四个参数就是这个世界的旋钮, 这也是科研仿真的日常:调参,看曲线,再调。(这四行的翻译是我们加的,方程和数值在书里。)

8. 原文地图¶

主题	原书章	原文位置
「难点在理解问题」点题	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:4`(搜「lies mainly in understanding」)
Navier–Stokes 与千禧年难题	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:9`(搜「Navier-Stokes」) · `text/29-ch13-13-numerical-mathematics.txt:19`(搜「millennial」)
风洞与原型的现实	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:34`(搜「wind tunnels」)
高斯曲线无解析积分、总面积为 1	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:53`(搜「Gaussian bell curve」)
黎曼和公式 13.2、两步实现	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:60`(搜「Riemann sum」) · `text/29-ch13-13-numerical-mathematics.txt:75`(搜「split the interval」)
中点切片技巧与 0.750	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:105`(搜「x_support」) · `text/29-ch13-13-numerical-mathematics.txt:129`(搜「0.750」)
「代码七行」框	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:131`(搜「Understanding the Problem」)
梯形法与 0.740	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:147`(搜「trapezoidal」) · `text/29-ch13-13-numerical-mathematics.txt:154`(搜「0.740」)
发散函数与误差上界	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:185`(搜「diverges」) · `text/29-ch13-13-numerical-mathematics.txt:190`(搜「upper bound」)
蒙特卡洛积分与区间变换	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:195`(搜「Monte Carlo」) · `text/29-ch13-13-numerical-mathematics.txt:222`(搜「multiplying」)
微分方程与谐振子方程	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:228`(搜「Differential Equations」) · `text/29-ch13-13-numerical-mathematics.txt:237`(搜「mẍ」)
解析解与「我们不感兴趣」	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:250`(搜「not interested」)
欧拉法公式与初始条件	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:265`(搜「Euler method」) · `text/29-ch13-13-numerical-mathematics.txt:266`(搜「initial conditions」)
第一、二步的手算	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:276`(搜「first time step」)
循环实现与画图对照	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:292`(搜「v = [v_0]」) · `text/29-ch13-13-numerical-mathematics.txt:310`(搜「x_theory」)
耦合方程与捕食者-猎物	13 Numerical Mathematics	`text/29-ch13-13-numerical-mathematics.txt:320`(搜「Coupled」) · `text/29-ch13-13-numerical-mathematics.txt:329`(搜「Ḃ = aB」) · `text/29-ch13-13-numerical-mathematics.txt:341`(搜「B = [10]」)

出处:「13 Numerical Mathematics」第 3–7 段(text/29-ch13-13-numerical-mathematics.txt:8,搜「no exact solution」)。 ↩
出处:「13 Numerical Mathematics」第 9–24 段(text/29-ch13-13-numerical-mathematics.txt:10,搜「basis of aerodynamics」;第 19 段,搜「million dollars」)。 ↩
出处:「13 Numerical Mathematics」第 26–34 段(text/29-ch13-13-numerical-mathematics.txt:34,搜「wind tunnels」)。 ↩
出处:「13 Numerical Mathematics」第 45–56 段(text/29-ch13-13-numerical-mathematics.txt:53,搜「cannot be integrated analytically」)。 ↩
出处:「13 Numerical Mathematics」第 58–129 段(text/29-ch13-13-numerical-mathematics.txt:75,搜「split the interval」;中点技巧见第 104–109 段;0.750 见第 129 段)。 ↩↩↩↩↩↩↩
出处:「13 Numerical Mathematics」第 131–141 段(加框「Understanding the Problem」,text/29-ch13-13-numerical-mathematics.txt:133,搜「quite normal」;「never write code」见第 139 段,搜「haven’t understood」)。 ↩
出处:「13 Numerical Mathematics」第 143–176 段(text/29-ch13-13-numerical-mathematics.txt:147,搜「trapezoidal rule」;权重结构见第 149 段;0.740 见第 153 段)。 ↩↩↩↩
出处:「13 Numerical Mathematics」第 177–194 段(text/29-ch13-13-numerical-mathematics.txt:185,搜「diverges」;误差上界公式 13.4 见第 190 段)。 ↩↩
出处:「13 Numerical Mathematics」第 195–227 段(text/29-ch13-13-numerical-mathematics.txt:210,搜「number_of_supporting_points」;区间变换与乘宽见第 218–226 段)。 ↩↩↩
出处:「13 Numerical Mathematics」第 228–234 段(text/29-ch13-13-numerical-mathematics.txt:232,搜「aerodynamics of aircraft」)。 ↩
出处:「13 Numerical Mathematics」第 235–251 段(text/29-ch13-13-numerical-mathematics.txt:237,搜「mẍ = −kx」;解析解与「不感兴趣」见第 244–250 段)。 ↩↩
出处:「13 Numerical Mathematics」第 253–290 段(text/29-ch13-13-numerical-mathematics.txt:265,搜「Euler method」;初始条件见第 269–274 段;第一、二步手算见第 276–284 段)。 ↩↩
出处:「13 Numerical Mathematics」第 291–318 段(text/29-ch13-13-numerical-mathematics.txt:291,搜「repeat the calculation in a loop」;解析解对照见第 305–317 段,搜「x_theory」)。 ↩↩↩↩↩
出处:「13 Numerical Mathematics」第 320–350 段(text/29-ch13-13-numerical-mathematics.txt:322,搜「predators and prey」;方程见第 329 段;代码见第 338–350 段)。 ↩↩↩↩↩↩