7.2 科学研究领域：从数值实验到复现科学

7.2 科学研究领域：从数值实验到可复现科学

本节摘要：科学计算把 Matlab 用作数值实验的通用底座，从求解微分方程到蒙特卡洛模拟。本节以偏微分方程数值解与随机模拟两类典型计算为主线，讲求解器使用，并给出可复现科学的三要素实践。

实验室的第二张桌子

理论推导有黑板，实验验证有台子，中间还有一张"数值实验桌"——在这张桌上，假设可以被快速计算检验，参数可以被系统扫描。Matlab 在科研中的位置就是这张桌子：交互环境让"改一个参数重算"的成本趋近于零，这正是假设检验需要的节奏。

典型计算一：微分方程数值解

ODE 求解器家族按刚度分工，非刚性问题用 ode45，刚性问题换 ode15s：

% 经典范德波振子
f = @(t,y) [y(2); (1 - y(1)^2)*y(2) - y(1)];
[t, y] = ode45(f, [0 20], [2; 0]);
plot(t, y(:,1)); xlabel('时间'); ylabel('位移');
​

求解器自动控制步长，使用者只管提供右端函数。刚性的判断有经验可循：系统里快慢动态差几个数量级（化学动力学、电路），ode45 会慢得无法忍受，换 ode15s 立竿见影——选错求解器不是错答案，是等不到答案。

PDE 则用有限元工作流：图形化界面画几何、设材料与边界条件、网格剖分、求解、后处理一气呵成，热传导与结构力学问题不必手写离散化。

完整案例：一次参数扫描研究

背景：非线性课程要求研究范德波振子的非线性强度参数如何改变振动形态。操作按"单点验证—扫描—汇总"三段推进。先在默认参数下确认解的形态无误，再扫描：

mus = [0.1 0.5 1 3 8];           % 非线性强度扫描表
feat = zeros(numel(mus), 2);     % 每档记录周期与幅值
for i = 1:numel(mus)
    mu = mus(i);
    f = @(t,y) [y(2); mu*(1 - y(1)^2)*y(2) - y(1)];
    [t, y] = ode45(f, [0 100], [2; 0]);
    % 峰值间隔的中位数作为"等效周期"
    pk = findpeaks(y(:,1));
    feat(i,:) = [median(diff(pk(:,1))), max(y(:,1))];
end
disp(table(mus.', feat(:,1), feat(:,2), ...
    'VariableNames', {'mu','周期','幅值'}))
​

结果解读：弱非线性（mu 等于 0.1）时波形近正弦、周期接近线性振子；mu 到 8 时出现典型的弛豫振荡——缓慢蓄能加快速释放，周期大幅拉长。一张参数—周期曲线就把"非线性如何重塑动力学"讲清楚，这比推导更早地给出直觉，也常是推导方向的来源。变式：把 ode45 换成 ode15s 对比耗时，mu 大时系统刚性渐显，两个求解器的效率差会直观显现——参数扫描顺手变成一次求解器选型实验。

ODE 求解器速查与常见事故

高频事故三则。一是右端函数里混入了绘图或打印，每步都执行一遍，仿真慢如蜗牛，右端必须纯净。二是容差设得过松导致轨迹漂移，关键物理量（能量、守恒量）不守恒时要收紧相对容差并检查右端是否写错符号。三是把代数环或事件漏掉：碰撞、开关类问题应使用事件机制在精确时刻中断积分，靠缩小步长硬碰只会得到含糊的过冲。科研代码的通病则是只调求解器不验网格——PDE 结果对网格的依赖没做收敛性检查之前，任何精巧的后处理都建在沙滩上。

科研侧还有个值得写进习惯的小纪律：结果与生成分离。分析脚本算出的每个数字、每张图，都应当能从"原始数据加脚本加种子"完整重生成，中间结果尽量不落盘。落盘的中间结果一旦被后续脚本悄悄依赖，"原始数据"就悄悄换成了"上一次运行的结果"，可复现链条由此断裂。需要缓存加速时，把缓存文件标记清楚并在脚本里可开关，是兼顾效率与纯净的折中。论文返修时能一键重跑全部图表的人，和面对审稿人质疑只能逐个截图翻找的人，体验过的是两种科研人生。

数值实验的记录还有个轻量格式可循：每个实验目录四个文件——入口脚本、参数记录、原始数据引用、结果图表，加一个两行的说明文件写"做了什么、看到了什么"。这套结构十分钟即可建好，却让半年后的你能在一分钟内找到任何一张图的出处。科研代码的组织成本远低于其回报，这是数值实验桌纪律的最后一块拼图。

典型计算二：蒙特卡洛模拟

rng(2026, 'twister');        % 固定随机流：可复现第一要素
n = 1e6;
est = zeros(1, 20);
for k = 1:20                  % 看估计量的波动
    p = rand(n, 2);
    est(k) = 4 * mean(sum(p.^2, 2) <= 1);
end
std(est)                      % 收敛速度按根号 n 下降
​

蒙特卡洛的误差按样本数的平方根倒数下降——精度提高一位数，样本要多两个数量级。知道收敛阶，才知道该加样本还是该改方法（方差缩减技术是后话）。

可复现三要素

"论文里的图重不出来"是学术界的陈年痼疾，Matlab 侧的解法是三件事：

• 版本：论文注明 Matlab 版本与工具箱清单（ver 的输出存档）；
• 随机种子：所有随机模拟开头 rng 固定种子，附在补充材料里；
• 数据快照：save 保存产生每张图的原始变量与脚本，一个图一个快照目录。

图 可复现数值实验的组成

💡 关键直觉：数值实验和物理实验一样需要"实验记录本"。快照目录就是记录本的一页：脚本、数据、种子、版本缺一不可，缺页的实验等于没做过。

本节要点回顾

• 数值实验桌：交互环境把"改参数重算"成本压到零，匹配假设检验的节奏；
• ODE 按刚性选求解器：ode45 通用、ode15s 治刚性；
• 蒙特卡洛收敛按平方根倒数，加样本前先想方差缩减；
• 可复现三要素：版本清单、随机种子、数据快照；
• 一图一档，把实验记录本的纪律带进代码。