翻开《核工程导论》（Lamarsh）或者《核反应堆分析》（Duderstadt & Hamilton），几乎每个核工学生都会被那种优雅的物理之美打动。中子输运方程的推导行云流水，扩散理论、临界分析、各种精巧的数学模型一字排开，你在纸面上推演得酣畅淋漓，仿佛已经掌握了反应堆的一切。

可一旦打开VS Code、PyCharm或者Jupyter Notebook，想要真正动手搭一个哪怕最简化的反应堆模拟器，你撞上的是一堵教科书极少提及的墙——巨大的工程鸿沟。

这个鸿沟并不在物理层面。方程还是那些方程，原理还是那些原理。真正的问题在于：你学到的那些连续、光滑、用微分写成的数学，计算机根本不认识。计算机只认识数字，而你面对的是一堆需要翻译成数值算法的微分方程、边界条件和反应率参数。翻译的过程，就是下面这张路线图：

连续物理过程 → 数值离散化 → 线性代数 → 算法设计 → 能跑的Python代码。大部分经典教材只走到第一步，科学软件的灵魂藏在后四步里。

第一步：把微分方程离散化，让计算机“看得见”

拿反应堆物理里最基础的方程来说——一维稳态中子扩散方程，带一个外中子源：

D·d²φ/dx² - Σₐφ + S = 0

这里面D是扩散系数，φ是中子通量，Σₐ是宏观吸收截面，S是外源强度。方程确实很漂亮，但微分符号d/dx对计算机来说形同天书。计算机只会做加减乘除，它理解不了导数。

解决的办法就是用有限差分近似，把连续的空间切成一个个等间距的节点。二阶导数的近似形式就变成了：

d²φ/dx² ≈ (φ_{i+1} - 2φ_i + φ_{i-1}) / Δx²

把这个近似塞回原方程，微分方程瞬间降级为一组代数方程：

D × (φ_{i+1} - 2φ_i + φ_{i-1}) / Δx² - Σₐφ_i + S_i = 0

到这里，反应堆不再是一个光滑的连续体，而是一排分别写着各自通量值的离散格子。你在做的事情，已经没有解析推导的优雅，而是实打实的“把几何掰碎”。

第二步：把代数关系堆成线性方程组

接下来要做的，是把上面那个节点方程重新排列，把含φ_{i-1}、φ_i、φ_{i+1}的项各自归拢。整理完以后，每个内部节点都遵循同一个模子：

(D/Δx²)·φ_{i-1} + ( -2D/Δx² - Σₐ )·φ_i + (D/Δx²)·φ_{i+1} = -S_i

这一下，物理问题彻底变了样：我们不再需要求解微分方程，而是要解一个巨大的线性方程组 AΦ = B。矩阵A是一个三对角矩阵，每行的非零元素只出现在对角线相邻的位置上，物理意义是相邻节点之间的中子泄漏和本节点的吸收损失；未知向量Φ就是每个节点上的中子通量；右端项B装着外源信息。

反应堆物理，走到这儿，已经悄然变成了数值线性代数。

第三步：用代码把矩阵“喂”给求解器

一旦模型被打包成矩阵形式，编程这件事反而变得异乎寻常地直接。下面是极简版的Python实现，用的是NumPy的线性方程求解器，几乎就是照着上面那套逻辑一行行码出来：

import numpy as npdef solve_1d_reactor_flux(core_length, num_nodes, D, Sigma_a, source_strength):    dx = core_length / (num_nodes - 1)    x = np.linspace(0, core_length, num_nodes)        A = np.zeros((num_nodes, num_nodes))    B = np.zeros(num_nodes)        leakage = D / dx**2    center = -(2*leakage + Sigma_a)        for i in range(1, num_nodes - 1):        A[i, i-1] = leakage        A[i, i]   = center        A[i, i+1] = leakage        B[i] = -source_strength        # 真空边界条件：两端通量固定为零    A[0, 0] = 1    A[-1, -1] = 1    B[0] = 0