用Python模拟三个经典科学实验：蒙特卡洛、随机游走与单摆数值积分

小时候在科学课上老师往一杯清水里轻轻放下一枚回形针水面竟然像一层薄薄的膜一样托住了金属把几滴牛奶滴进盘子再蘸一点洗洁精颜色就会迅速四散开。这些“哇”的一瞬间背后往往藏着表面张力、分子运动、能量转换这些硬核原理。后来开始写代码才发现很多看起来“只能靠肉眼观察”的小实验其实完全可以用 Python 还原出来而且模拟的过程能让人把原理看得更清楚、更精确。这篇文章就用三个十分经典的趣味实验来展开用随机撒点估算圆周率、用随机游走模拟布朗运动、用数值积分模拟单摆摆动。你不需要有高深的物理背景也不需要很长的代码经验只要能运行 Python就可以跟着一步步把现象“造”出来。读完你不仅能得到三份可复现的代码更能理解“现象背后的数学和物理”以及程序中容易踩的坑。1. 背景与核心概念1.1 什么是趣味科学小实验趣味科学小实验通常是指那些用简单道具就能完成、现象直观、结果让人惊讶的小型实验。它最大的特点是“门槛低但原理不浅”。比如一个单摆只需要一根绳子和一个重物但它的运动方程背后是二阶常微分方程再比如往正方形里随机撒豆子看起来是小朋友的游戏但用到的却是蒙特卡洛方法和概率论。这类实验非常适合作为编程练手项目因为它们的物理规律明确、数学表达清晰同时又不像大型工程项目那样需要复杂架构。用 Python 实现时你可以把注意力集中在“如何把物理规律写成代码”“如何让结果可视化”“如何处理数值误差”这些问题上。1.2 为什么用 Python 做趣味实验用 Python 做科学小实验有几点天然优势。首先是语法简洁。一个蒙特卡洛实验的核心逻辑可能只需要十几行代码新手可以快速看到效果不会因为语言细节而失去耐心。其次是生态完善。numpy 提供高效的数组运算matplotlib 可以快速绘图scipy 还能解微分方程。这三个库基本覆盖了从数据生成、数值计算到可视化展示的完整链路。第三是容易可视化。科学实验最讲究“观察现象”而 matplotlib 能轻松画出随机点的分布、粒子的运动轨迹、角度随时间的变化曲线。现象可视化之后原理的理解难度会降低很多。1.3 本文实验总览本文挑选了三个有代表性的实验实验名称核心知识点代码难度蒙特卡洛法估算圆周率随机抽样、几何概率、大数定律入门布朗运动模拟随机游走、扩散定律、统计平均初级单摆数值模拟常微分方程、欧拉法、Runge-Kutta 法中级这三个实验从简单到复杂恰好覆盖了“随机模拟”“统计可视化”“数值积分”这几个在工程和科研中非常常见的编程场景。学会它们你不仅能做这三个实验还能迁移到许多更复杂的模拟任务中。2. 环境准备与版本说明2.1 运行环境本文代码以 Python 3.8 为基准如果你的环境是 Python 3.10 或 3.11 也完全没有问题。代码中主要使用以下第三方库numpy用于数组计算、随机数生成、累加求和。matplotlib用于绘制散点图、折线图展示实验结果。操作系统不限Windows、macOS、Linux 都可以运行。本文示例以常见的命令行环境为例如果你使用的是 PyCharm、VS Code 等 IDE同样适用。版本不需要刻意追求最新因为这三个实验只用到了 numpy 和 matplotlib 最基础、最稳定的功能。重点在于理解实验思路而不是追新版本。2.2 安装依赖库建议先创建一个独立目录并在该目录下创建虚拟环境避免污染系统全局环境。mkdir python_lab cd python_lab python -m venv venv激活虚拟环境Windowsvenv\Scripts\activatemacOS / Linuxsource venv/bin/activate然后安装依赖pip install numpy matplotlib安装完成后可以快速验证一下python -c import numpy, matplotlib; print(ok)如果输出ok说明环境已经准备好了。2.3 项目结构为了便于管理建议按实验编号创建三个文件python_lab/ ├── venv/ ├── exp1_monte_carlo_pi.py ├── exp2_brownian_motion.py └── exp3_pendulum.py每个文件都是独立的可以直接运行。后续如果要扩展可以在这个结构上继续添加新的实验文件。3. 实验一蒙特卡洛法估算圆周率3.1 实验现象先想象一个场景一张纸上画一个边长为 2 的正方形正方形内部画一个半径为 1 的圆。现在你闭着眼睛往纸上撒豆子统计落在圆内的豆子数量再统计总豆子数量。理论上圆面积和正方形面积之比是 π/4所以圆内豆子比例应该也接近 π/4。于是π 的估计值就可以写成π ≈ 4 × (圆内豆子数) / (总豆子数)这就是著名的蒙特卡洛估算 π 的方法。它并不需要真的去测量圆的周长只需要“随机撒点”和“数数”。3.2 原理分析蒙特卡洛方法的核心思想是用大量随机样本来估算一个确定性的数值。当样本数量足够多时频率会趋近于概率。在本实验中随机点(x, y)在正方形内均匀分布。判断点在圆内的条件是x² y² ≤ 1圆面积与正方形面积之比为圆面积 / 正方形面积 πr² / (2r)² π/4所以只要统计圆内点的比例再乘以 4就得到了 π 的近似值。这里有一个非常重要的隐含条件随机样本必须均匀分布。如果随机数生成器有偏或者采样范围不对估算结果就会出错。3.3 完整代码# 文件路径exp1_monte_carlo_pi.py import random import matplotlib.pyplot as plt def estimate_pi(num_points: int, seed: int None): 使用蒙特卡洛方法估算圆周率。 参数 num_points: 随机点数量 seed: 随机种子传入后结果可复现 返回 pi_estimate: π 的估算值 points_in: 圆内点坐标列表 points_out: 圆外点坐标列表 if seed is not None: random.seed(seed) inside 0 points_in [] points_out [] for _ in range(num_points): x random.uniform(-1, 1) y random.uniform(-1, 1) if x * x y * y 1: inside 1 points_in.append((x, y)) else: points_out.append((x, y)) pi_estimate 4.0 * inside / num_points return pi_estimate, points_in, points_out def plot_result(points_in, points_out, pi_estimate): 将圆内外的点绘制到图像中。 fig, ax plt.subplots(figsize(6, 6)) if points_in: xs, ys zip(*points_in) ax.scatter(xs, ys, s1, colorblue, label圆内点) if points_out: xs, ys zip(*points_out) ax.scatter(xs, ys, s1, colorred, label圆外点) circle plt.Circle((0, 0), 1, fillFalse, colorblack, linewidth1.5) ax.add_patch(circle) ax.set_aspect(equal) ax.set_xlim(-1.2, 1.2) ax.set_ylim(-1.2, 1.2) ax.legend() ax.set_title(f蒙特卡洛估算 π ≈ {pi_estimate:.4f}) plt.show() if __name__ __main__: N 10000 pi_value, inside_points, outside_points estimate_pi(N, seed42) print(f点数{N}估算 π 值{pi_value:.4f}) plot_result(inside_points, outside_points, pi_value)3.4 运行与结果分析运行代码python exp1_monte_carlo_pi.py你会看到一个正方形和圆的散点图蓝色点落在圆内红色点落在圆外。同时命令行会输出一个 π 的估算值。在seed42、点数10000的条件下结果通常会落在3.12到3.16之间而真实 π 是3.14159...。这说明即使只有一万个样本估算结果也已经比较接近了。你可以试着修改N的值当N1000时波动会更大有时候可能是3.08有时候可能是3.18。当N100000时结果会更稳定前两位小数大概率能稳定在3.14附近。这个现象背后的理论依据是“大数定律”样本量越大频率越接近概率。但也要注意蒙特卡洛方法给出的结果本身带有随机误差它的收敛速度通常是 O(1/√N)也就是说每提高一位小数精度所需点数大约要增加一百倍。4. 实验二布朗运动随机游走模拟4.1 实验现象在显微镜下观察悬浮在水中的花粉颗粒你会发现它们一直在做无规则运动一会儿向左一会儿向右方向完全不可预测。这个现象就是布朗运动它是由水分子从四面八方撞击花粉颗粒造成的。用 Python 模拟布朗运动最经典的方式是“随机游走”让一个粒子每一步都随机选择一个方向然后不断累积位置。虽然单条轨迹看起来毫无规律但如果你同时模拟大量粒子并统计它们距离起点的平均平方距离就会发现一个非常漂亮的统计规律。4.2 原理分析为了简化模拟我们可以把二维平面上的每一步看成x 方向随机选择1或-1y 方向随机选择1或-1。那么每走一步粒子位移平方的期望值是E[x² y²] E[x²] E[y²] 1 1 2如果粒子走了 t 步由于每一步相互独立均方位移Mean Squared Displacement, MSD满足MSD 2t也就是说尽管单条轨迹看起来完全随机但大量粒子的统计平均结果却呈现出“位移平方随时间线性增长”的规律。这正是爱因斯坦在 1905 年解释布朗运动时得到的核心结论也是扩散现象的微观基础。4.3 完整代码# 文件路径exp2_brownian_motion.py import numpy as np import matplotlib.pyplot as plt def simulate_particles(n_particles: int 200, steps: int 500, seed: int None): 模拟多个粒子的二维随机游走。 返回 positions: 形状为 (n_particles, steps1, 2) 的轨迹数组 msd: 每个时间步上所有粒子的均方位移形状为 (steps1,) time: 时间序列范围 [0, steps] if seed is not None: np.random.seed(seed) # 每步在 x 和 y 方向各自随机选择 1 或 -1 step_choices np.random.choice([-1.0, 1.0], size(n_particles, steps, 2)) # 沿着时间轴累加得到每个粒子的绝对位置 positions np.cumsum(step_choices, axis1) # 在最前面补一个原点表示初始位置 init np.zeros((n_particles, 1, 2)) positions np.concatenate([init, positions], axis1) # 均方位移每个时刻对所有粒子的 x²y² 求平均 msd np.mean(np.sum(positions ** 2, axis2), axis0) time np.arange(steps 1) return positions, msd, time def plot_brownian(positions, msd, time): 绘制前 10 条粒子的轨迹以及均方位移曲线。 fig, (ax1, ax2) plt.subplots(1, 2, figsize(12, 5)) # 左图部分粒子的运动轨迹 for i in range(10): ax1.plot(positions[i, :, 0], positions[i, :, 1], linewidth0.8) ax1.set_title(前 10 个粒子的随机游走轨迹) ax1.set_xlabel(x) ax1.set_ylabel(y) ax1.grid(alpha0.3) # 右图均方位移随时间变化 ax2.plot(time, msd, label模拟 MSD, linewidth2) ax2.plot(time, 2 * time, linestyle--, label理论值 MSD 2t, linewidth2) ax2.set_title(均方位移随时间变化) ax2.set_xlabel(步数 t) ax2.set_ylabel(MSD) ax2.legend() ax2.grid(alpha0.3) plt.tight_layout() plt.show() if __name__ __main__: positions, msd, time simulate_particles(n_particles500, steps400, seed2024) # 观察最后一步的模拟值与理论值 theoretical 2 * time[-1] print(f最后一步模拟 MSD{msd[-1]:.2f}) print(f理论 MSD{theoretical:.2f}) plot_brownian(positions, msd, time)4.4 运行与结果分析运行代码python exp2_brownian_motion.py你会看到左右两张图。左边是一团乱麻一样的随机轨迹每条轨迹都完全不同右边是一条近似直线的 MSD 曲线理论直线MSD 2t与模拟结果基本重合。这说明一个很重要的思想单个随机事件无法预测但大量随机事件汇总之后会体现出稳定的统计规律。这也是为什么“随机”并不等于“没有规律”。可以继续修改n_particles和steps来观察变化粒子数量太少时MSD 曲线会抖动明显和理论直线偏差较大。粒子数量足够多时比如 1000 个以上MSD 曲线会非常平滑。这里的模拟结果不仅适用于布朗运动还能推广到股票价格随机波动、动物觅食路径、群体扩散等许多场景。5. 实验三单摆数值模拟与积分误差5.1 实验现象单摆是中学物理里最常见的模型一根轻绳一端固定另一端挂一个小球轻轻推动小球它就会周期性地摆动。很多人会想当然地认为“单摆的运动就是正弦函数”但实际上这只在小角度近似下才成立。真实单摆的运动方程是非线性的当摆角较大时运动规律和正弦曲线有明显差别。更神奇的地方在于程序的数值求解过程同样的物理问题用不同的数值积分方法结果可能差很多。实验中我会用欧拉法和四阶 Runge-Kutta 法RK4分别求解同样的单摆方程对比角度的变化曲线以及系统能量的变化。你可以清楚地看到即使是看起来“很简单”的单摆数值方法选不好也会得到完全错误的结果。5.2 原理分析忽略空气阻力单摆的运动方程为d²θ/dt² -(g / L) * sin(θ)其中 θ 是摆角g 是重力加速度L 是摆长。为了写成程序我们引入角速度 ωdθ/dt ω dω/dt -(g / L) * sin(θ)欧拉法的更新方式非常直观θ_{n1} θ_n h * ω_n ω_{n1} ω_n - h * (g / L) * sin(θ_n)这种方法的优点是简单缺点是误差累积较快尤其在大摆角、大步长时系统的能量会逐渐漂移甚至出现“越摆越高”这种违背物理规律的现象。四阶 Runge-Kutta 方法的思路是用四个不同位置的斜率加权平均从而更准确地估计每一步的变化量。虽然公式多一些但精度远高于欧拉法。对于长时间模拟RK4 是更稳的选择。5.3 完整代码# 文件路径exp3_pendulum.py import numpy as np import matplotlib.pyplot as plt G 9.8 L 1.0 def derivatives(state): 计算单摆状态 (theta, omega) 的导数。 theta, omega state dtheta omega domega -(G / L) * np.sin(theta) return np.array([dtheta, domega]) def euler_step(state, dt): 欧拉法单步更新。 theta, omega state dtheta, domega derivatives(state) return np.array([theta dt * dtheta, omega dt * domega]) def rk4_step(state, dt): 四阶 Runge-Kutta 单步更新。 k1 derivatives(state) k2_state state 0.5 * dt * k1 k2 derivatives(k2_state) k3_state state 0.5 * dt * k2 k3 derivatives(k3_state) k4_state state dt * k3 k4 derivatives(k4_state) return state (dt / 6.0) * (k1 2 * k2 2 * k3 k4) def simulate(methodrk4, dt0.01, total_time20.0, theta01.0): 从头开始模拟单摆运动。 steps int(total_time / dt) times np.linspace(0, total_time, steps 1) state np.array([theta0, 0.0]) states np.zeros((steps 1, 2)) states[0] state for i in range(1, steps 1): if method euler: state euler_step(state, dt) else: state rk4_step(state, dt) states[i] state return times, states def mechanical_energy(state): 计算单摆机械能。 取小球质量 m1能量公式为 E 0.5 * L² * ω² g * L * (1 - cosθ) theta, omega state kinetic 0.5 * L * L * omega * omega potential G * L * (1 - np.cos(theta)) return kinetic potential if __name__ __main__: DT 0.02 TOTAL_TIME 30.0 THETA0 1.2 times_euler, states_euler simulate(euler, DT, TOTAL_TIME, THETA0) times_rk4, states_rk4 simulate(rk4, DT, TOTAL_TIME, THETA0) energy_euler [mechanical_energy(s) for s in states_euler] energy_rk4 [mechanical_energy(s) for s in states_rk4] # 小角度近似解析解仅用于参考对比 analytic_theta THETA0 * np.cos(np.sqrt(G / L) * times_rk4) fig, (ax1, ax2) plt.subplots(2, 1, figsize(12, 8)) # 上半部分角度随时间变化 ax1.plot(times_euler, states_euler[:, 0], label欧拉法, linewidth1.5) ax1.plot(times_rk4, states_rk4[:, 0], labelRK4 法, linewidth1.5) ax1.plot(times_rk4, analytic_theta, linestyle--, label小角度近似解析解, linewidth1.5) ax1.set_xlabel(时间 t/s) ax1.set_ylabel(摆角 θ/rad) ax1.legend() ax1.grid(alpha0.3) # 下半部分机械能随时间变化 ax2.plot(times_euler, energy_euler, label欧拉法能量, linewidth1.5) ax2.plot(times_rk4, energy_rk4, labelRK4 法能量, linewidth1.5) ax2.set_xlabel(时间 t/s) ax2.set_ylabel(机械能 E) ax2.legend() ax2.grid(alpha0.3) plt.tight_layout() plt.show() print(30 秒模拟结束后的欧拉法能量, energy_euler[-1]) print(30 秒模拟结束后的 RK4 法能量, energy_rk4[-1])5.4 运行与结果分析运行代码python exp3_pendulum.py你会看到上下两个子图。上方是摆角随时间的变化下方是机械能随时间的变化。观察结果通常有以下特征欧拉法在dt0.02、模拟 30 秒的条件下角度峰值会缓慢变化能量曲线明显上升或下降说明能量不再守恒。RK4 法角度曲线非常平滑能量基本维持在初始值附近系统看起来更接近真实物理。小角度近似解析解与 RK4 法在大摆角θ01.2时会有明显差异这是非线性效应导致的。这个实验最大的启示是程序运行“成功”不等于结果“正确”。对于数值模拟必须选择合适的方法和步长并用能量守恒等物理约束来验证结果是否可信。这一点在做工程仿真、游戏物理、机器人控制时尤其重要。6. 常见问题与排查思路6.1 matplotlib 中文乱码在 Matplotlib 绘图时如果标题、图例中出现中文很可能显示成方框或乱码。这是因为 Matplotlib 默认字体不支持中文。解决办法是显式指定中文字体比如在代码开头加上import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, Arial Unicode MS] plt.rcParams[axes.unicode_minus] False不同操作系统可用的字体不完全一样需要根据本机已安装的中文字体调整。如果你不想处理字体也可以把图中的中文标签全部换成英文。6.2 随机数导致结果不稳定蒙特卡洛实验每次运行结果都不同这是正常的并不是程序出错。要保证结果可复现只需要在生成随机数之前固定随机种子。例如在实验一中传入seed42在实验二中传入seed2024再运行几次结果就会保持一致。6.3 单摆能量不守恒如果你用欧拉法模拟单摆并计算系统的机械能会发现能量要么越来越大要么越来越小。这是因为欧拉法属于低阶显式格式数值耗散问题严重。解决办法使用 RK4 等更高阶的方法减小步长dt还可以改用 symplectic integrator 等专门针对物理系统设计的积分器。在工程应用中判断一个仿真程序是否可靠不能只看界面是否正常更要用能量、动量等守恒量做“合理性检查”。6.4 其他常见问题速查问题现象常见原因解决思路ModuleNotFoundError: No module named numpy没有安装依赖执行pip install numpy matplotlib图形窗口一闪而过脚本运行结束后窗口被关闭在代码末尾加plt.show()或在脚本末尾加input()暂停单摆角度越来越大欧拉法步长太大减小dt或改用 RK4散点图太密集看不清点数太多抽样绘制部分点或者把点的大小调小程序能运行但图像空白变量未正确传入绘图函数检查点列表是否为空打印长度确认7. 最佳实践与工程建议7.1 使用随机种子保证可复现只要实验里出现随机数就要在代码里考虑可复现性。最直接的办法是提供seed参数例如本文实验一和实验二中的写法。这样在调试、对拍、写文档时都能稳定复现同一组结果。在团队协作中可复现性比“偶尔得到一个好结果”重要得多。你写下的每一个随机种子都是在告诉未来的自己“这个结果是基于什么条件得到的”。7.2 把计算逻辑和可视化分离观察本文三个实验的代码可以发现我都把“计算”和“绘图”拆成了不同函数。这是一种很低成本但收益很高的习惯estimate_pi只负责计算和返回数据plot_result只负责把数据画出来。这样拆分之后你可以非常方便地换一种可视化方式或者把计算结果导出成 CSV甚至接入其他系统而不需要改动核心算法。7.3 数值实验的性能优化实验一用 for 循环逐点判断优点是直观缺点是点数变大时速度变慢。实际处理大规模随机采样时更推荐用 numpy 实现向量化版本。以下是实验一的向量化改写思路import numpy as np def estimate_pi_vectorized(num_points: int, seed: int None): if seed is not None: np.random.seed(seed) points np.random.uniform(-1, 1, size(num_points, 2)) distances points[:, 0] ** 2 points[:, 1] ** 2 inside int(np.sum(distances 1)) pi_estimate 4.0 * inside / num_points return pi_estimate同样的逻辑用 for 循环写十万个点会明显变慢而 numpy 版本可以轻松处理百万甚至千万个点。这就是向量化带来的性能红利。7.4 实验记录与验证做数值实验时建议每次运行都记录以下内容随机种子关键参数比如摆长、步长、总时长、采样点数输出结果比如估算的 π 值和 MSD运行环境比如 Python 版本和依赖库版本。这样做的目的是让实验可以追溯。一旦发现结果异常你能很快定位是参数问题、代码问题还是环境问题。对于正式项目而言这种记录习惯几乎是必须的。8. 总结与下一步学习到这里我们完成了三个趣味小实验的完整模拟。你不仅仅看到了现象还掌握了背后的原理用蒙特卡洛方法估算 π理解了“大量随机样本可以逼近确定值”的思想用随机游走模拟布朗运动观察到“单粒子随机群体统计却呈现规律”的扩散定律用欧拉法和 RK4 模拟单摆理解了数值积分精度和能量校验的重要性。下一步你可以尝试把这些实验继续做深把蒙特卡洛实验改成可视化“估算误差随点数变化”的收敛曲线把布朗运动换成三维随机游走观察均方位移与维度的关系把单摆扩展成双摆体验混沌现象感受初始条件的敏感性试试用scipy.integrate.solve_ivp替代手写求解器对比不同方法的误差。如果本篇对你有帮助建议收藏备用尤其是三个实验的完整代码和排错表格在实际写代码的时候可以随时翻出来参考。动手跑一遍再试着改一改参数你会发现“代码里的物理”比想象中更有趣。