数学建模与Python:高空坠落的硬币传说

from numpy import sqrtfrom pint import UnitRegistry# 创建单位注册表units = UnitRegistry()# 定义基本单位meter = units.metersecond = units.second

# 定义参数a = 9.8 * meter / second**2# 重力加速度h = 381 * meter              # 帝国大厦高度# 计算落地时间t = sqrt(2 * h / a)print(f"落地时间: {t:.2f}")# 计算落地速度v = a * tprint(f"落地速度: {v:.2f}")

mile = units.milehour = units.hourv_mph = v.to(mile/hour)print(f"时速: {v_mph:.2f}")

v_terminal = 29 * meter / secondv_terminal_mph = v_terminal.to(mile/hour)print(f"实际落地时速: {v_terminal_mph:.2f}")

import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationimport matplotlib.patches as patches# ==========================================# 物理模型参数# ==========================================H = 381.0# 帝国大厦高度 (m)g = 9.8# 重力加速度 (m/s^2)v_term = 29.0# 终端速度 (m/s)# 时间设置dt = 0.1t_max = 15.0# 模拟时长steps = int(t_max / dt)time = np.linspace(0, t_max, steps)# 1. 真空模型 (Vacuum Model)y_vac = H - 0.5 * g * time**2v_vac = g * time# 落地修正mask_vac = y_vac < 0y_vac[mask_vac] = 0v_vac[mask_vac] = 0# 2. 空气阻力模型 (Air Resistance Model)# 简化模拟：先加速到终端速度，然后匀速# 达到终端速度的时间t_term = v_term / gy_air = np.zeros_like(time)v_air = np.zeros_like(time)for i, t in enumerate(time):if t < t_term:# 加速阶段        v_air[i] = g * t        y_air[i] = H - 0.5 * g * t**2else:# 匀速阶段        v_air[i] = v_term# y = y_at_term - v_term * (t - t_term)        y_at_term = H - 0.5 * g * t_term**2        y_air[i] = y_at_term - v_term * (t - t_term)# 落地修正if y_air[i] < 0:        y_air[i] = 0        v_air[i] = 0# ==========================================# 炫酷可视化设置# ==========================================plt.style.use('dark_background')fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6), gridspec_kw={'width_ratios': [1, 1.5]})# --- 左图：大楼全景视角 (Empire State Building View) ---ax1.set_xlim(-10, 20)ax1.set_ylim(0, 420)ax1.set_title("Empire State Drop Test", color='white', fontsize=14, pad=20)ax1.set_ylabel("Height (m)")# 绘制大楼轮廓building = patches.Rectangle((-5, 0), 10, 381, linewidth=2, edgecolor='#3498db', facecolor='#2c3e50', alpha=0.5)ax1.add_patch(building)# 地面ground = ax1.axhline(0, color='#2ecc71', lw=3)ax1.text(-8, 5, 'Ground', color='#2ecc71', fontsize=10)# 硬币对象coin_vac, = ax1.plot([], [], 'o', color='#e74c3c', markersize=8, label='Vacuum Model') # 红色coin_air, = ax1.plot([], [], 'o', color='#f1c40f', markersize=8, label='Real Model')   # 黄色# 轨迹线trail_vac, = ax1.plot([], [], '--', color='#e74c3c', alpha=0.5, lw=1)trail_air, = ax1.plot([], [], '-', color='#f1c40f', alpha=0.5, lw=1)ax1.legend(loc='upper right', frameon=False)# --- 右图：速度仪表盘 & 冲击力模拟 (Dashboard View) ---ax2.set_xlim(0, t_max)ax2.set_ylim(0, 100)ax2.set_title("Velocity & Impact Analysis", color='white', fontsize=14, pad=20)ax2.set_xlabel("Time (s)")ax2.set_ylabel("Velocity (m/s)")ax2.grid(True, linestyle='--', alpha=0.3)# 速度曲线line_v_vac, = ax2.plot([], [], '--', color='#e74c3c', lw=2, label='Vacuum Velocity')line_v_air, = ax2.plot([], [], '-', color='#f1c40f', lw=2, label='Real Velocity')# 终端速度参考线ax2.axhline(v_term, color='#f1c40f', linestyle=':', alpha=0.6)ax2.text(0.5, v_term + 2, f'Terminal Velocity: {v_term} m/s', color='#f1c40f', fontsize=8)# 实时数据文本text_h = ax2.text(0.05, 0.9, '', transform=ax2.transAxes, color='white', fontsize=10)text_v = ax2.text(0.05, 0.85, '', transform=ax2.transAxes, color='white', fontsize=10)# 落地累积效果 (Impact Accumulation)impact_scatter = ax1.scatter([], [], color='#ecf0f1', s=20, alpha=0.6, marker='x')impact_points = []# 镜头切换逻辑# 0-5s: 全景; 5-10s: 追踪下落; 10s+: 落地特写defupdate(frame):    current_time = time[frame]# 1. 更新位置和速度    h_v = y_vac[frame]    h_a = y_air[frame]    vel_v = v_vac[frame]    vel_a = v_air[frame]# 更新硬币    coin_vac.set_data([0], [h_v])    coin_air.set_data([5], [h_a])# 更新轨迹    trail_vac.set_data([0]*frame, y_vac[:frame])    trail_air.set_data([5]*frame, y_air[:frame])# 更新速度曲线    line_v_vac.set_data(time[:frame], v_vac[:frame])    line_v_air.set_data(time[:frame], v_air[:frame])# 更新文本    text_h.set_text(f"Height: Vac={h_v:.1f}m | Real={h_a:.1f}m")    text_v.set_text(f"Speed:  Vac={vel_v:.1f}m/s | Real={vel_a:.1f}m/s")# 2. 镜头控制 (Camera Control)if h_a > 100:# 高空俯瞰        ax1.set_ylim(h_a - 200, h_a + 50)elif h_a > 0:# 接近地面        ax1.set_ylim(0, 200)else:# 落地后保持地面视角        ax1.set_ylim(0, 50)# 3. 落地堆积效果if frame % 5 == 0: # 模拟多次实验的累积             impact_points.append([np.random.normal(5, 1), 0])if len(impact_points) > 0:                 arr = np.array(impact_points)                 impact_scatter.set_offsets(arr)return coin_vac, coin_air, trail_vac, trail_air, line_v_vac, line_v_air, text_h, text_v, impact_scatter# 生成动画ani = FuncAnimation(fig, update, frames=steps, interval=30, blit=True)ani.save('01_falling_penny_cool.gif', writer='pillow', fps=30)plt.close(fig)

2. 动态演示效果

通过这个动画，我们可以清晰地观察到：

1. 初始阶段：两者几乎同步下落。
2. 加速阶段：真空模型（红色）持续加速，迅速超越现实模型（黄色）。
3. 终端速度：现实模型在大约 3 秒后达到终端速度（29 m/s），之后匀速下落。
4. 最终结果：真空中的硬币像子弹一样砸向地面，而现实中的硬币则缓缓飘落（相对而言），并没有传说中那么可怕。

总结

通过这个案例，我们学到了建模的核心思想：

1. 所有模型都是不准确的：模型一忽略了空气阻力，它是不准确的。模型二假设硬币瞬间达到终端速度，也是简化的，严格来说也是有误差的。
2. 但有些模型是有用的：模型二虽然简化，但它给出的结论（硬币砸不死人）与《流言终结者》的实验结果一致。这就足够解决我们的问题了。

这就是科学建模的魅力：我们不需要复刻整个宇宙，只需要捕捉到关键的相关规律。