Гравитационная задача N тел
Пример показывает бесконечную 2D-анимацию системы из 64 взаимно
притягивающихся тел. Тело с индексом 0 намного массивнее остальных, но оно не
закреплено: получает ускорение от других тел и движется вместе со всей системой.
Модель
Для каждого тела хранятся масса, координаты, скорость и ускорение. На каждом
шаге сначала обнуляются ускорения, затем один раз обрабатывается каждая пара
i, j, где i < j:
dx := x[j] - x[i]
dy := y[j] - y[i]
r2 := dx * dx + dy * dy + SOFTENING2
invR3 := 1.0 / (sqrt(r2) * r2)
ax[i] := ax[i] + G * mass[j] * dx * invR3
ay[i] := ay[i] + G * mass[j] * dy * invR3
ax[j] := ax[j] - G * mass[i] * dx * invR3
ay[j] := ay[j] - G * mass[i] * dy * invR3
Оба тела пары обновляются одновременно равными и противоположными силами. Это полноценное взаимное N-body взаимодействие, а не набор независимых орбит в фиксированном центральном поле.
SOFTENING2 ограничивает ускорение при близком сближении тел. Без softening
дискретный интегратор потребовал бы очень малого шага времени.
Интегрирование velocity Verlet
Начальные ускорения вычисляются до входа в цикл. На каждом шаге velocity Verlet сначала обновляет координаты с ускорением в начале шага:
x[i] := x[i] + vx[i] * DT + 0.5 * ax[i] * DT * DT
y[i] := y[i] + vy[i] * DT + 0.5 * ay[i] * DT * DT
После этого ускорения пересчитываются по новым координатам, а скорость корректируется средним ускорением в начале и конце шага:
vx[i] := vx[i] + 0.5 * (ax[i] + nextAx[i]) * DT
vy[i] := vy[i] + 0.5 * (ay[i] + nextAy[i]) * DT
Velocity Verlet существенно лучше сохраняет энергию орбитальной системы, чем явный или semi-implicit Euler при том же шаге времени.
Начальная скорость массивного тела компенсирует суммарный импульс остальных, поэтому центр масс системы не получает искусственный начальный дрейф.
Визуализация
У каждого тела сохраняются последние TAIL положений. Тело 0 рисуется
крупнее, остальные получают разные оттенки. Вызов новый лист завершает кадр и
подготавливает следующий чёрный лист.
Полный исходный код находится в
examples/painter/orbits.kum.