积分具有二次阻力的抛体
在落地事件处停止，与真空参考比较，并用阻力功核对机械能损失。
目标与运行方法
本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。
学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。
环境与参数
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from scipy.integrate import solve_ivp
from IPython.display import display as _display
def display(value):
if isinstance(value,pd.DataFrame):
table=value.to_html(max_rows=12,escape=True)
style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
return _display(value)
RECIPE_ID='ph03-projectile'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
"axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6 Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads Core data mode: embedded snapshot; no network requests
读取内嵌快照
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
"""Read embedded data from this notebook, including a renamed copy."""
for path in sorted(Path.cwd().glob("*.ipynb")):
try:
notebook = json.loads(path.read_text(encoding="utf-8"))
metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
blob = metadata.get("snapshots", {}).get(snapshot_key)
if metadata.get("recipe_id") != recipe_id or not blob:
continue
if blob.get("sha256") != expected_sha256:
continue
raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
if hashlib.sha256(raw).hexdigest() != expected_sha256:
raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
return json.loads(raw)
except (OSError, json.JSONDecodeError):
continue
raise FileNotFoundError(
"Save/open the complete .ipynb in the active working folder. "
"The embedded snapshot could not be located for " + recipe_id)
SNAPSHOT_HASHES={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf', 'constants': '2533802a1a546c612d06ba2e931e4f32a418be854c75a5e9c327ff57d5b56126'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['models', 'constants']
方法函数
def projectile_with_drag(mass,diameter,drag_coefficient,air_density,speed,angle_deg,gravity=9.80665,height_m=1.5):
"""Constant-density quadratic drag; terminate at a descending ground crossing."""
if mass<=0 or diameter<=0 or drag_coefficient<0 or air_density<0 or speed<=0 or gravity<=0 or height_m<0:
raise ValueError('Invalid projectile parameters')
area=np.pi*(diameter/2)**2;factor=.5*air_density*drag_coefficient*area/mass
angle=np.deg2rad(angle_deg);initial=[0,height_m,speed*np.cos(angle),speed*np.sin(angle)]
def dynamics(t,state):
x,z,vx,vz=state;velocity=np.hypot(vx,vz)
return [vx,vz,-factor*velocity*vx,-gravity-factor*velocity*vz]
def ground(t,state):return state[1]
ground.direction=-1;ground.terminal=True
solved=solve_ivp(dynamics,[0,60],initial,events=ground,dense_output=True,rtol=1e-9,atol=1e-11,max_step=.1)
if not solved.success or len(solved.t_events[0])!=1:raise RuntimeError('Descending ground event not found')
flight=float(solved.t_events[0][0]);t=np.linspace(0,flight,501);state=solved.sol(t)
velocity=np.hypot(state[2],state[3]);energy=.5*mass*velocity**2+mass*gravity*state[1]
drag_power=-.5*air_density*drag_coefficient*area*velocity**3
return dict(t=t,state=state,flight_s=flight,range_m=float(state[0,-1]),energy_J=energy,drag_power_W=drag_power)
print("Method ready: projectile_with_drag")
Method ready: projectile_with_drag
数据与模型边界
单摆、振子、抛体和温度轨迹是明确声明的数值模型，不是实验室测量。使用的 g0 为约定标准重力加速度。GW150914 使用真实冻结应变和已发表事件时间；滤波曲线不是独立的显著性或源参数分析。黑体曲线假定理想辐射体。
分步分析
1. 设定力模型及落地事件
使用恒定空气密度、阻力系数及球形迎风截面，省略风、旋转及速度相关系数。事件在下降的 z=0 过零处停止，撞击后的积分不再属于飞行轨迹。
p=snapshot("models")["projectile"];gravity=snapshot("constants")["constants"]["g_0"]["value"]
arguments=dict(mass=p["mass_kg"],diameter=p["diameter_m"],air_density=p["air_density_kg_per_m3"],speed=p["speed_m_per_s"],angle_deg=p["angle_deg"],gravity=gravity,height_m=1.5)
drag=projectile_with_drag(**arguments,drag_coefficient=p["drag_coefficient"])
vacuum=projectile_with_drag(**arguments,drag_coefficient=0)
fig,ax=plt.subplots(figsize=(8,4));ax.plot(drag["state"][0],drag["state"][1],label="Quadratic drag")
ax.plot(vacuum["state"][0],vacuum["state"][1],"--",label="Vacuum")
ax.set(xlabel="Horizontal displacement (m)",ylabel="Height (m)",title="Declared flight model until first descending impact");ax.legend();plt.tight_layout();plt.show()
2. 检查速度与飞行时间
速率为 sqrt(vx²+vz²)，单位 m/s。两条曲线在各自撞击时间结束，不能假定持续时间相等而强行对齐。此例中的阻力同时影响轨迹和飞行时间。
fig,ax=plt.subplots(figsize=(8,4))
for name,solution in [("Drag",drag),("Vacuum",vacuum)]:ax.plot(solution["t"],np.hypot(solution["state"][2],solution["state"][3]),label=name)
ax.set(xlabel="Flight time (s)",ylabel="Speed (m/s)",title="Speed before impact");ax.legend();plt.tight_layout();plt.show()
display(pd.DataFrame({"model":["drag","vacuum"],"flight_s":[drag["flight_s"],vacuum["flight_s"]],"range_m":[drag["range_m"],vacuum["range_m"]]}))
| model | flight_s | range_m | |
|---|---|---|---|
| 0 | drag | 3.889601 | 60.470905 |
| 1 | vacuum | 4.395881 | 93.250713 |
3. 用阻力功核对机械能损失
机械能为 ½mv²+mgz，单位焦耳。阻力功率为负，单位瓦。对功率做时间积分，应在离散精度内恢复能量变化。这是声明模型中的守恒审计，不是阻力实测。
energy_change=drag["energy_J"]-drag["energy_J"][0]
widths=np.diff(drag["t"]);power=drag["drag_power_W"]
work=np.r_[0,np.cumsum(widths*(power[:-1]+power[1:])/2)]
fig,ax=plt.subplots(figsize=(8,4));ax.plot(drag["t"],energy_change,label="Mechanical-energy change")
ax.plot(drag["t"],work,"--",label="Integrated drag work")
ax.set(xlabel="Flight time (s)",ylabel="Energy change / work (J)",title="Energy balance of the dissipative model");ax.legend();plt.tight_layout();plt.show()
4. 检查真空极限并导出
零阻力参考具有解析飞行时间和射程，应包含非零发射高度。解释阻力差异之前，检查最终落地高度和能量平衡。
angle=np.deg2rad(arguments["angle_deg"]);vx=arguments["speed"]*np.cos(angle);vz=arguments["speed"]*np.sin(angle)
expected_time=(vz+np.sqrt(vz*vz+2*gravity*arguments["height_m"]))/gravity
assert np.isclose(vacuum["flight_s"],expected_time,rtol=1e-8)
assert np.isclose(vacuum["range_m"],vx*expected_time,rtol=1e-8)
assert abs(drag["state"][1,-1])<1e-8 and np.all(np.diff(drag["energy_J"])<=1e-5)
assert np.isclose(energy_change[-1],work[-1],rtol=.002)
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"time_s":drag["t"],"x_m":drag["state"][0],"z_m":drag["state"][1],"vx_m_per_s":drag["state"][2],"vz_m_per_s":drag["state"][3],"mechanical_energy_J":drag["energy_J"],"drag_power_W":power}).to_csv(OUTPUT_DIR/"drag-flight.csv",index=False)
print("Ground event, analytic vacuum and work checks passed; exports:",OUTPUT_DIR)
Ground event, analytic vacuum and work checks passed; exports: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads/results/ph03-projectile
自主练习与参考答案
- 说明阻力功率的符号。
- 为何真空参考要包含发射高度？
参考解释。 阻力反向于速度，其与速度的内积非正。非零高度会改变飞行时间与射程。
先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。
可选的 AI 编程练习
提示词 1
仅用 NumPy、SciPy、Matplotlib 在同一阻力模型下比较三十、四十五、六十度。
提示词 2
使用同样的包改变恒定空气密度，导出轨迹及能量/功误差。
人工验收。 在下降落地事件停止，零阻力核对及有量纲能量平衡通过。
检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。
可选联网扩展
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://physics.nist.gov/cuu/Constants/',timeout=20)
response.raise_for_status()
print("Source:",response.url)
print(response.text[:700])
except requests.RequestException as error:
print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged
来源、快照与下一步
NIST CODATA physical constants — NIST public reference data; credit NIST/CODATA; retrieved 2026-10-09. SHA-256
77fb90e66c40db3e6eb16630bc9c88e4c7c8beddbe5e71be406f2f26e3f67e67. Preserve official ASCII table; extract seven explicitly named SI constants and their units.IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256
a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256
4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.
下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。