速度与数值积分

已知每个位置的风速，如何画出粒子的下一位置？如果速度始终为 10 m/s，60 秒的位移是 600 m。速度随位置和时间变化时，可以把大段时间切成小步，每步重新查询速度，再累积位移。这就是数值积分（numerical integration）的基本想法。

一阶导数与积分相互连接

dx/dt = v(x,t) 表示位置随时间的局部变化率等于速度。Euler 法（Euler method）使用本步起点的速度，假设它在本小步中保持不变：

$$\frac{d\mathbf{x}}{dt}=\mathbf{v}(\mathbf{x},t),{\mathbf{x}}_{n+1}={\mathbf{x}}_n+\mathbf{v}({\mathbf{x}}_n,t_n)\mathrm{\Delta }t$$

符号 / 英文	含义	单位
x / position	平面位置，或带明确定义的坐标	本段为 m
v / velocity	速度向量	m/s
Δt / time step	一次推进的模拟时长	s
dx/dt / time derivative	位置相对于时间的变化率	m/s
ω / angular velocity	绕中心旋转的角速度	rad/s

如果一步迈得太远，就会跳过速度变化明显的区域。屏幕帧率高并不自动保证模拟准确：时间倍率很大时，每帧对应的模拟步长仍可能很长。

手算：绕圆运动为何慢慢跑偏

在平面设 v(x,y)=(-ωy,ωx)，其中 ω=0.01 rad/s。它处处沿圆的切线方向，理想轨迹保持离原点的距离不变。粒子从 (1000,0) m 出发，取 Δt=10 s：

起点速度 = (0,10) m/s
Euler下一点 = (1000,100) m
离中心距离 = sqrt(1000²+100²) ≈ 1004.99 m

一次就向外偏了约 5 m。原因是使用起点切线代替了整段弧线；减小步长能减小这种离散误差，但不会让有限步长的 Euler 变为精确圆周运动。

中点法（midpoint method，一种二阶 Runge–Kutta / RK2）先预测半步位置，再用中点的速度走完整一步：

type XY = readonly [number, number]
const velocity = ([x, y]: XY): XY => [-0.01*y, 0.01*x]
function midpointStep(p: XY, dt: number): XY {
  const a = velocity(p)
  const mid: XY = [p[0] + a[0]*dt/2, p[1] + a[1]*dt/2]
  const b = velocity(mid)
  return [p[0] + b[0]*dt, p[1] + b[1]*dt]
}
midpointStep([1000, 0], 10) // [995,100]，半径约1000.0125 m

它改善了本例误差，但仍非精确解。若速度场随时间变化，第一次查询用 t，第二次应使用 t+Δt/2。更高阶方法仍依赖正确的坐标、单位、采样和步长。

U/V 怎样更新经纬度

先采用球形地球近似，半径 R，纬度 φ、经度 λ，U/V 分别为局部东向和北向速度。沿纬圈和经线的小位移对应：

$$\mathrm{\Delta }\lambda \approx \frac{U\mathrm{\Delta }t}{R\cos \phi },\mathrm{\Delta }\phi \approx \frac{V\mathrm{\Delta }t}{R}$$

右边单位为弧度。北纬 60° 向东 10 m/s，60 s 后经度约增加 0.01079°。这不是给经度加 600，而是先将米制位移转为角位移。极区 cosφ 接近零，经度参数退化；本简化公式适合远离极区的小步长区域演示，精密全球推进应选择适合的椭球或三维方法。

本书风粒子实验实际使用区域球面近似和 Euler 推进。这里的 RK2 是独立数学例子，说明改进方向，并不意味着实验提供了算法切换。

自检：时间倍率与粒子数

原来每一步为 1 s，将模拟速度调为 100 倍时，能否直接改成每步 100 s？粒子数从 1024 增至 4096，速度应如何变化？

答案与理由

可以产生某种动画，但误差可能显著增加。更稳妥的方式是增加若干子步，并限制每步穿过的网格距离；具体阈值取决于场变化与精度要求。粒子数只改变采样密度和计算成本，不应改变物理速度。还应分别记录真实经过的时间、展示倍率和实际完成的模拟时间，避免卡顿时误读结果。

下一步：精度、误差与数值检查。回到主课：GPU 风粒子。