省流

运动的描述是求解动力学过程的基础。拉格朗日量和哈密顿量的求解十分依赖坐标的选取。一套好的坐标可能隐含了约束，从而降低求解的难度。

固定坐标系下的质点运动学

矢量力学通过科学抽象，将研究对象简化为刚体乃至质点。我们通过矢量来描述质点的诸多物理量：

选定参考系的原点，质点的位置由径向矢量 $\vec{r}=r\hat{\tau}$ 表示；
质点的运动通过径向矢量的变化——位移来表示：$\Delta \vec{r}=\vec{r}_2-\vec{r}_1$；
质点的平均速度由位移对时间的导数来表示：$\bar{\vec{v}}=\frac{\Delta \vec{r}}{\Delta t}$；
质点的平均加速度由速度对时间的导数来表示：$\bar{\vec{a}}=\frac{\Delta \vec{v}}{\Delta t}$。

很多时候，我们更关心质点的瞬时速度和瞬时加速度：

$\begin{aligned} \vec{v}&=\lim_{\Delta t \to 0} \frac{\Delta \vec{r}}{\Delta t}=\frac{d\vec{r}}{dt}=\dot{r}\hat{\tau}+r\dot{\hat{\tau}}\\ \vec{a}&=\lim_{\Delta t \to 0} \frac{\Delta \vec{v}}{\Delta t}=\frac{d\vec{v}}{dt}=\ddot{r}\hat{\tau}+2\dot{r}\dot{\hat{\tau}}+r\ddot{\hat{\tau}}\\ \end{aligned}$

代入$\dot{\hat{\tau}}=\dot{\theta}\hat{n},\dot{\hat{n}}=-\dot{\theta}\hat{\tau}$，其中$\hat{n}$为极坐标中的横向单位矢量：

$\begin{aligned} \vec{v}&=\dot{r}\hat{\tau}+r\dot{\theta}\hat{n}\\ \vec{a}&=\ddot{r}\hat{\tau}+2\dot{r}\dot{\theta}\hat{n}+r\ddot{\theta}\hat{n}-r\dot{\theta}^2\hat{\tau}\\ \end{aligned}$

如果不把$\vec{v}$按径向横向分解，而是按轨道本身的切向与法向分解，可以得到：
$\begin{aligned} \vec{a}&=\dot{v}\hat{e}_t+v\dot{\theta}\hat{e}_n\\ &=\dot{v}\hat{e}_t+v\frac{d\theta}{ds}\frac{ds}{dt}\hat{e}_n\\ &=\dot{v}\hat{e}_t+\frac{v^2}{\rho}\hat{e}_n\\ \end{aligned}$
其中$\hat{e}_t$为单位切向矢量，$\hat{e}_n$为单位法向矢量（指向曲率中心，也称主法向），$\rho$为曲率半径。

一般我们会在具体的坐标系下进行计算，具体的坐标系会影响基矢的选取。

直角坐标系

在直角坐标系中：

$\begin{aligned} \vec{r}&=x\hat{i}+y\hat{j}+z\hat{k}\\ \vec{v}&=\dot{x}\hat{i}+\dot{y}\hat{j}+\dot{z}\hat{k}\\ \vec{a}&=\ddot{x}\hat{i}+\ddot{y}\hat{j}+\ddot{z}\hat{k}\\ \end{aligned}$

极坐标系

在极坐标系中：

$\begin{aligned} \vec{r}&=r\hat{\tau}\\ \vec{v}&=\dot{r}\hat{\tau}+r\dot{\theta}\hat{n}\\ \vec{a}&=\ddot{r}\hat{\tau}+2\dot{r}\dot{\theta}\hat{n}+r\ddot{\theta}\hat{n}-r\dot{\theta}^2\hat{\tau}\\ \end{aligned}$

其中，速度分为径向速度$\dot{r}$和切向速度$r\dot{\theta}$，加速度分为径向加速度$\ddot{r}-r\dot{\theta}^2$和切向加速度$r\ddot{\theta}+2\dot{r}\dot{\theta}$。其中径向加速度包括：

径向速度的大小变化：$\ddot{r}$；
向心加速度（横向速度的方向变化）：$-r\dot{\theta}^2$。

切向加速度包括：

横向速度的大小变化：$r\ddot{\theta}$；
径向速度的方向变化：$2\dot{r}\dot{\theta}$。

柱坐标系

在柱坐标系中：

$\begin{aligned} \vec{r}&=r\hat{\tau}+z\hat{k}\\ \vec{v}&=\dot{r}\hat{\tau}+r\dot{\theta}\hat{n}+\dot{z}\hat{k}\\ \vec{a}&=\ddot{r}\hat{\tau}+2\dot{r}\dot{\theta}\hat{n}+r\ddot{\theta}\hat{n}-r\dot{\theta}^2\hat{\tau}+\ddot{z}\hat{k}\\ \end{aligned}$

柱坐标系的基矢只有$\hat{\tau},\hat{n}$随角度$\theta$转动，$\hat{k}$不变。在直角坐标系中$\hat{\tau}=(\cos\theta,\sin\theta,0),\ \hat{n}=(-\sin\theta,\cos\theta,0)$，对时间求导即得：
$\dot{\hat{\tau}}=\dot{\theta}\hat{n},\quad \dot{\hat{n}}=-\dot{\theta}\hat{\tau},\quad \dot{\hat{k}}=0$

球坐标系

在球坐标系中：

$\begin{aligned} \vec{r}&=r\hat{\tau}\\ \vec{v}&=\dot{r}\hat{\tau}+r\dot{\theta}\hat{\theta}+r\dot{\phi}\sin\theta\hat{\phi}\\ \vec{a}&=(\ddot{r}-r\dot{\theta}^2-r\dot{\phi}^2\sin^2\theta)\hat{\tau}\\&+(2\dot{r}\dot{\theta}+r\ddot{\theta}-r\dot{\phi}^2\sin\theta\cos\theta)\hat{\theta}\\&+(r\ddot{\phi}\sin\theta+2\dot{r}\dot{\phi}\sin\theta+2r\dot{\theta}\dot{\phi}\cos\theta)\hat{\phi}\\ \end{aligned}$

推导该表达式，最关键的是推导基矢的导数。一方面可以通过几何图像直接得出；另一方面可以在直角坐标系下写出各基矢，对时间求导后再反代回球坐标系。两种方法都得到：
$\begin{aligned} \dot{\hat{\tau}}&=\dot{\theta}\hat{\theta}+\dot{\phi}\sin\theta\,\hat{\phi}\\ \dot{\hat{\theta}}&=-\dot{\theta}\hat{\tau}+\dot{\phi}\cos\theta\,\hat{\phi}\\ \dot{\hat{\phi}}&=-\dot{\phi}\sin\theta\,\hat{\tau}-\dot{\phi}\cos\theta\,\hat{\theta} \end{aligned}$
以$\dot{\hat{\tau}}$为例演示直角坐标法：先在直角坐标下写出基矢$\hat{\tau}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$，对时间求全导数（$\theta,\phi$均为$t$的函数，用链式法则）：
$\dot{\hat{\tau}}=\left(\dot\theta\cos\theta\cos\phi-\dot\phi\sin\theta\sin\phi,\ \dot\theta\cos\theta\sin\phi+\dot\phi\sin\theta\cos\phi,\ -\dot\theta\sin\theta\right)$
按$\dot\theta,\dot\phi$拆项并提取公因子，再对照直角坐标表达式识别基矢（反解）：
$\dot{\hat{\tau}}=\dot\theta(\cos\theta\cos\phi,\cos\theta\sin\phi,-\sin\theta)+\dot\phi\sin\theta(-\sin\phi,\cos\phi,0)=\dot\theta\hat{\theta}+\dot\phi\sin\theta\hat{\phi}$
更系统的方法是利用基矢的正交归一性，把$\dot{\hat{\tau}}$按基矢展开，系数即点积：$\dot{\hat{\tau}}\cdot\hat{\tau}=\frac12\frac{d}{dt}|\hat{\tau}|^2=0$（单位矢量垂直于自己的导数）、$\dot{\hat{\tau}}\cdot\hat{\theta}=\dot\theta$、$\dot{\hat{\tau}}\cdot\hat{\phi}=\dot\phi\sin\theta$，结果相同。$\hat{\theta},\hat{\phi}$的导数用同一流程（分别对$(\cos\theta\cos\phi,\cos\theta\sin\phi,-\sin\theta)$和$(-\sin\phi,\cos\phi,0)$求全导数再投影回基矢）。将$\vec{r}=r\hat{\tau}$对时间连续求导两次并代入基矢导数，即得速度与加速度的各分量。

另一种方法是通过拉格朗日方程：
$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}\right)=\frac{\partial L}{\partial q},\quad L=\frac12m\left(\dot{r}^2+r^2\dot{\theta}^2+r^2\sin^2\theta\,\dot{\phi}^2\right)-V$
来推出：
$\begin{cases}\dfrac{d}{dt}\left(\dfrac{\partial L}{\partial \dot{r}}\right)&=m\ddot{r}\\ \dfrac{\partial L}{\partial {r}}&=mr(\dot{\theta}^2+\sin^2\theta\,\dot{\phi}^2)-\dfrac{\partial V}{\partial r}\end{cases}\Rightarrow a_r=\ddot{r}-r(\dot{\theta}^2+\sin^2\theta\,\dot{\phi}^2)$
与直接求导的结果一致。其中，拉格朗日量中的动能项来自速度的平方：在正交坐标系下速度平方的交叉项为零，动能退化为各项平方和（对角元），即$T=\frac12m(\dot{r}^2+r^2\dot{\theta}^2+r^2\sin^2\theta\,\dot{\phi}^2)$，这也是拉格朗日力学中动能总是广义速度二次型的原因。

相对运动

描述运动的时候总会选择一个方便的参考系，但是到了求解具体问题的时候，根据问题选择的参考系显然更简便，这就需要不同参考系之间的互相转换。

已知质点相对于运动参考系$S’$的位置$\vec{r}’$、速度$\vec{v}’$、加速度$\vec{a}’$。另一个的参考系$S$相对于运动参考系的位置、速度、加速度分别为$\vec{R}$、$\vec{V}$、$\vec{A}$，当我们站在这个参考系内，我们认为这个参考系是固定的，这样会方便我们对问题的建模。

任意平动参考系

在新参考系的表示为：

$\begin{cases} \vec{r}(t)&=\vec{r}'+\vec{R}(t)\\ \vec{v}(t)&=\vec{v}'+\vec{V}(t)\\ \vec{a}(t)&=\vec{a}'+\vec{A}(t) \end{cases}$

在不涉及转动（即基矢变化）的新坐标系中，我们发现所有物理量都可以做简单平移来转化得到新的坐标系下的物理量。

旋转参考系

取时间$t=0$的时候两个坐标系基矢重合，则该时间零点瞬间，新坐标系下的物理量表达为：

$\begin{cases} \vec{r}&=\vec{r}'\\ \vec{v}&=\vec{v}'+\vec{\omega}\times\vec{r}\\ \vec{a}&=\vec{a}'+\vec{\beta}\times\vec{r}+2\vec{\omega}\times\vec{v}'+\vec{\omega}\times(\vec{\omega}\times\vec{r})\\ \end{cases}$

这里不区分两个坐标系下的位置矢量，因为在$t=0$时刻两个坐标系的原点重合，所以位置矢量相同。

这里给出证明：

$\begin{aligned} \vec{a}(t)&=\frac{d\vec{v}(t)}{dt}\\ &=\frac{d}{dt}(\vec{v}'+\vec{\omega}\times\vec{r})\\ &=\frac{d}{dt}\sum(v_i\vec{i}'+\epsilon_{ijk}\omega_i r_j\vec{k}')\\ &=\sum[(a_i'\vec{i}'+v_i'\vec{\omega}\times\vec{i}')+(\epsilon_{ijk}\beta_i r_j\vec{k}'+\epsilon_{ijk}\omega_i v_j'\vec{k}'+\epsilon_{ijk}\omega_i r_j\vec{\omega}\times\vec{k}')]\\ &=\vec{a}'+\vec{\beta}\times\vec{r}+2\vec{\omega}\times\vec{v}'+\vec{\omega}\times(\vec{\omega}\times\vec{r})\\ \end{aligned}$

这里也有更优美的证明。利用以下结论：
$(\frac{d}{dt}\vec{A})_{Still}=(\frac{d}{dt}\vec{A})_{Rotating}+\vec{\omega}\times\vec{A}$
可以得到：
$\begin{aligned} \vec{v}&=\frac{d\vec{r}}{dt}\\ &=\vec{v}'+\vec{\omega}\times\vec{r} \end{aligned}$ $\begin{aligned} \vec{\beta}&=\frac{d\vec{\omega}}{dt}\\ &=\frac{d}{dt}\vec{\omega}'+\vec{\omega}\times\vec{\omega}\\ &=\vec{\beta}' \end{aligned}$
代入下式：
$\begin{aligned} \vec{a}&=\frac{d\vec{v}}{dt}\\ &=\frac{d}{dt}(\vec{v}'+\vec{\omega}\times\vec{r})+\vec{\omega}\times(\vec{v}'+\vec{\omega}\times\vec{r})\\ &=\vec{a}'+\vec{\beta}'\times\vec{r}+2\vec{\omega}\times\vec{v}'+\vec{\omega}\times(\vec{\omega}\times\vec{r})\\ \end{aligned}$

三项牵连加速度，分别是：

$\vec{\beta}\times\vec{r}$：由于参考系的旋转加速度引起的加速度（角加速度）；
$2\vec{\omega}\times\vec{v}’$：由于相对参考系的运动引起的加速度（科里奥利加速度）；
$\vec{\omega}\times(\vec{\omega}\times\vec{r})$：由于参考系的旋转引起的加速度（向心加速度）。

实际上应用的时候，是因为我们想计算非惯性系中的加速度才使用相对运动理论，所以不知道的一般是$\vec{a}’$，这说明我们要把上式移项变成：

$\vec{a}'=\vec{a}-\vec{\beta}\times\vec{r}-2\vec{\omega}\times\vec{v}'-\vec{\omega}\times(\vec{\omega}\times\vec{r})$

例子：一个质点做匀速圆周运动（半径$r$、速率$v=\omega r$），如果在相对静止参考系（即同样转速的旋转参考系）下观察，则$\vec{v}’=0$。由：
$\vec{a}'=\vec{a}-2\vec{\omega}\times\vec{v}'-\vec{\omega}\times(\vec{\omega}\times\vec{r})$
其中$\vec{a}=-\omega^2 r\,\hat{r}=-\omega v\,\hat{r}$（$\hat{r}$为径向单位矢量），$-2\vec{\omega}\times\vec{v}’=0$，$-\vec{\omega}\times(\vec{\omega}\times\vec{r})=+\omega^2 r\,\hat{r}=+\omega v\,\hat{r}$，故$\vec{a}’=0$，能观察到物体静止。如果是一个质点本身静止（$\vec{a}=0$），但是在旋转参考系中观察，则$\vec{v}’=-\vec{\omega}\times\vec{r}$，此时：
$-2\vec{\omega}\times\vec{v}'=-2\omega^2 r\,\hat{r}=-2\omega v\,\hat{r},\quad -\vec{\omega}\times(\vec{\omega}\times\vec{r})=+\omega^2 r\,\hat{r}=+\omega v\,\hat{r}$
故$\vec{a}’=-2\omega v\,\hat{r}+\omega v\,\hat{r}=-\omega v\,\hat{r}=-\omega^2 r\,\hat{r}$，在旋转参考系中观察到物体正常旋转，有向心加速度。

例子2：在赤道上抛出一个质点，取东为$x$正向、北为$y$正向、上为$z$正向，地球自转角速度$\vec{\omega}=\omega\hat{y}$。上升阶段（$v_z>0$）科里奥利加速度$a_x=-2\omega v_z<0$，使物体获得向西的速度分量；下降阶段（$v_z<0$）科里奥利加速度反向，恰好抵消此前速度分量的继续增加，最终物体全过程向西偏移。可以很简单地推出理论公式：
$v_z(t)=v_0-gt,\quad a_x(t)=-2\omega v_z(t),\quad v_x(t)=\int a_x(t)dt=-2\omega(v_0t-\frac{1}{2}gt^2),\quad x(t)=\int v_x(t)dt=-\omega(v_0t^2-\frac{1}{3}gt^3)$
代入$t=2v_0/g$，可以得到落点处向西的位移为：
$x(t)=-\omega\left(v_0\left(\frac{2v_0}{g}\right)^2-\frac{1}{3}g\left(\frac{2v_0}{g}\right)^3\right)=-\frac{4\omega v_0^3}{g^2}+\frac{8\omega v_0^3}{3g^2}=-\frac{4\omega v_0^3}{3g^2}$

科里奥利力

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 常量
omega = 7.292e-5  # 地球自转角速度 (rad/s)
g = 10          # 重力加速度 (m/s²)
v0 = 50.0         # 初始上抛速度 (m/s)
h0 = 0.0          # 初始高度 (m)
lat = 0.0         # 赤道纬度 (rad)

# 运动方程：dy/dt = f(t, y), y = [x, y, z, vx, vy, vz]
def motion_eq(t, y):
    x, y, z, vx, vy, vz = y
    
    # 科里奥利加速度 (a_c = -2 ω × v)
    a_cx = -2 * omega * (-np.sin(lat) * vy + np.cos(lat) * vz)
    a_cy = -2 * omega * np.sin(lat) * vx
    a_cz = -2 * omega * np.cos(lat) * vx
    
    # 总加速度 (重力 + 科里奥利)
    ax = a_cx
    ay = a_cy
    az = -g + a_cz
    
    return [vx, vy, vz, ax, ay, az]

# 初始条件 [x, y, z, vx, vy, vz]
y0 = [0.0, 0.0, h0, 0.0, 0.0, v0]

# 模拟时间 (上抛至落地共需 2*v0/g = 10 秒)
t_span = (0, 10)  # 10 秒足够上抛下落
t_eval = np.linspace(0, 10, 1000)

# 解微分方程
sol = solve_ivp(motion_eq, t_span, y0, t_eval=t_eval, rtol=1e-6)

# 提取轨迹
x, y, z = sol.y[0], sol.y[1], sol.y[2]

# 找到落地点 (z 首次回落到 ≤ 0 的时刻)
landing_idx = np.argmax(z <= 0)
if z[landing_idx] > 0:
    landing_idx = len(z) - 1
landing_x, landing_y = x[landing_idx], y[landing_idx]
theory = -4/3*omega*v0**3/g**2

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(x, z, label='Trajectory (x-z plane)')
plt.scatter(landing_x, 0, color='red', label=f'Landing: x={landing_x:.4f} m')
plt.axhline(0, color='black', linestyle='--', alpha=0.5)
plt.axvline(theory, color='red', linestyle='--', alpha=0.5,label='Theoretical Landing')
plt.xlabel('East-West Displacement (m)')
plt.ylabel('Height (m)')
plt.title('Projectile Motion with Coriolis Effect (Equator)')
plt.legend()
plt.grid()
plt.show()

# 输出结果
print(f"Initial position: (0, 0, {h0})")
print(f"Landing position: ({landing_x:.6f}, {landing_y:.6f}, 0)")
print(f"Westward displacement: {landing_x:.6f} m")