小时候在科学课上，老师往一杯清水里轻轻放下一枚回形针，水面竟然像一层薄薄的膜一样托住了金属；把几滴牛奶滴进盘子，再蘸一点洗洁精，颜色就会迅速四散开。这些“哇”的一瞬间，背后往往藏着表面张力、分子运动、能量转换这些硬核原理。后来开始写代码，才发现很多看起来“只能靠肉眼观察”的小实验，其实完全可以用 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激活虚拟环境：
- Windows：
venv\Scripts\activate- macOS / Linux：
source 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, s=1, color="blue", label="圆内点") if points_out: xs, ys = zip(*points_out) ax.scatter(xs, ys, s=1, color="red", label="圆外点") circle = plt.Circle((0, 0), 1, fill=False, color="black", linewidth=1.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, seed=42) print(f"点数：{N}，估算 π 值：{pi_value:.4f}") plot_result(inside_points, outside_points, pi_value)3.4 运行与结果分析
运行代码：
python exp1_monte_carlo_pi.py你会看到一个正方形和圆的散点图，蓝色点落在圆内，红色点落在圆外。同时命令行会输出一个 π 的估算值。
在seed=42、点数10000的条件下，结果通常会落在3.12到3.16之间，而真实 π 是3.14159...。这说明，即使只有一万个样本，估算结果也已经比较接近了。
你可以试着修改N的值：
- 当
N=1000时，波动会更大，有时候可能是3.08，有时候可能是3.18。 - 当
N=100000时，结果会更稳定，前两位小数大概率能稳定在3.14附近。
这个现象背后的理论依据是“大数定律”：样本量越大，频率越接近概率。但也要注意，蒙特卡洛方法给出的结果本身带有随机误差，它的收敛速度通常是 O(1/√N)，也就是说，每提高一位小数精度，所需点数大约要增加一百倍。
4. 实验二：布朗运动随机游走模拟
4.1 实验现象
在显微镜下观察悬浮在水中的花粉颗粒，你会发现它们一直在做无规则运动，一会儿向左，一会儿向右，方向完全不可预测。这个现象就是布朗运动，它是由水分子从四面八方撞击花粉颗粒造成的。
用 Python 模拟布朗运动，最经典的方式是“随机游走”：让一个粒子每一步都随机选择一个方向，然后不断累积位置。虽然单条轨迹看起来毫无规律，但如果你同时模拟大量粒子，并统计它们距离起点的平均平方距离，就会发现一个非常漂亮的统计规律。
4.2 原理分析
为了简化模拟，我们可以把二维平面上的每一步看成：
- x 方向随机选择
+1或-1； - y 方向随机选择
+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, steps+1, 2) 的轨迹数组 msd: 每个时间步上所有粒子的均方位移，形状为 (steps+1,) 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, axis=1) # 在最前面补一个原点，表示初始位置 init = np.zeros((n_particles, 1, 2)) positions = np.concatenate([init, positions], axis=1) # 均方位移：每个时刻，对所有粒子的 x²+y² 求平均 msd = np.mean(np.sum(positions ** 2, axis=2), axis=0) 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], linewidth=0.8) ax1.set_title("前 10 个粒子的随机游走轨迹") ax1.set_xlabel("x") ax1.set_ylabel("y") ax1.grid(alpha=0.3) # 右图：均方位移随时间变化 ax2.plot(time, msd, label="模拟 MSD", linewidth=2) ax2.plot(time, 2 * time, linestyle="--", label="理论值 MSD = 2t", linewidth=2) ax2.set_title("均方位移随时间变化") ax2.set_xlabel("步数 t") ax2.set_ylabel("MSD") ax2.legend() ax2.grid(alpha=0.3) plt.tight_layout() plt.show() if __name__ == "__main__": positions, msd, time = simulate_particles(n_particles=500, steps=400, seed=2024) # 观察最后一步的模拟值与理论值 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(θ)欧拉法的更新方式非常直观：
θ_{n+1} = θ_n + h * ω_n ω_{n+1} = ω_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(method="rk4", dt=0.01, total_time=20.0, theta0=1.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): """ 计算单摆机械能。 取小球质量 m=1，能量公式为： 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="欧拉法", linewidth=1.5) ax1.plot(times_rk4, states_rk4[:, 0], label="RK4 法", linewidth=1.5) ax1.plot(times_rk4, analytic_theta, linestyle="--", label="小角度近似解析解", linewidth=1.5) ax1.set_xlabel("时间 t/s") ax1.set_ylabel("摆角 θ/rad") ax1.legend() ax1.grid(alpha=0.3) # 下半部分：机械能随时间变化 ax2.plot(times_euler, energy_euler, label="欧拉法能量", linewidth=1.5) ax2.plot(times_rk4, energy_rk4, label="RK4 法能量", linewidth=1.5) ax2.set_xlabel("时间 t/s") ax2.set_ylabel("机械能 E") ax2.legend() ax2.grid(alpha=0.3) plt.tight_layout() plt.show() print("30 秒模拟结束后的欧拉法能量：", energy_euler[-1]) print("30 秒模拟结束后的 RK4 法能量：", energy_rk4[-1])5.4 运行与结果分析
运行代码：
python exp3_pendulum.py你会看到上下两个子图。上方是摆角随时间的变化，下方是机械能随时间的变化。
观察结果通常有以下特征：
- 欧拉法：在
dt=0.02、模拟 30 秒的条件下，角度峰值会缓慢变化，能量曲线明显上升或下降，说明能量不再守恒。 - RK4 法：角度曲线非常平滑，能量基本维持在初始值附近，系统看起来更接近真实物理。
- 小角度近似解析解与 RK4 法在大摆角
θ0=1.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 随机数导致结果不稳定
蒙特卡洛实验每次运行结果都不同，这是正常的，并不是程序出错。要保证结果可复现，只需要在生成随机数之前固定随机种子。
例如在实验一中传入seed=42，在实验二中传入seed=2024，再运行几次，结果就会保持一致。
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替代手写求解器，对比不同方法的误差。
如果本篇对你有帮助，建议收藏备用，尤其是三个实验的完整代码和排错表格，在实际写代码的时候可以随时翻出来参考。动手跑一遍，再试着改一改参数，你会发现“代码里的物理”比想象中更有趣。
