A. Noullez, D. Fanelli, E. Aurell による Heap-based algorithm for one-dimensional particle systems という論文を読んだ際のメモです.

はじめに

1次元 $N$ 体系の数値シミュレーションについて考えます. 1次元重力多体系を想定し, 個々の粒子の軌道は他の粒子と衝突しなければ解析的に求まると仮定します. しかし衝突 (または交差) が起こると粒子の配置が入れ替わりますから, 常にどこで衝突が起こるかを追跡し続ける必要があります. ナイーブにはこれには $O ( N )$ の時間がかかるように思えますが, この論文は $O ( \ln N )$ でよいと主張しています.

ヒープ

本記事において $R$ 分木が完全 (complete) であるとは, すべての枝の深さが $D$ または $D-1$ であり, 深さ $D$ の枝が可能な限り「左」に寄せられているもののことを言います. そのうえで, ヒープ とは, 各ノードにある値 (PartialOrd を要求します) が付随した完全 $R$ 分木であり, 常に親ノードの値はそのどの小ノードの値よりも小さい (または等しい) もののことを言います¹. 完全 $R$ 分木では, あるノードのインデックスから親ノードと子ノードのインデックスが決まっているため, ノード間の関係を保持する必要がなく, 実装が非常に簡単です. このあたり自分で実装し始めてから既に良さそうな Rust 実装 dary_heap があることに気づきましたが, この手のものはフルスクラッチで書かないと身につかないですし最後まで自分で実装しました.

ヒープの構築はナイーブには $O ( N \ln N )$ の時間がかかりますが, 効率的なアルゴリズムが存在し $O ( N )$ で可能なようです². いずれにせよ, 本記事のアルゴリズムではヒープは最初に一回構築すればもう一から作り直すことはしないので問題はないでしょう.

アルゴリズム

$N$ 体シミュレーションの最初に初期条件に従って粒子群を生成します. これは長さ $N$ の Vec<_> に確保されるでしょう. 次いで, 粒子 $j$ と粒子 $j+1$ が交差する時刻 $\tau_j$ のリストを生成します (これは長さ $N - 1$ です). 時刻 $\tau_j$ と粒子 id の組のリスト Vec<(usize, f64)> などとするのが良いと思います. この交差時刻のリストからヒープを構築します.

ヒープの定義から, そのルート要素は最も早い衝突時刻とそれがどの粒子で起こるかを教えてくれます. そこで, この粒子の組について, その位置と速度を衝突時刻まで進め (このとき当然両者の空間座標は一致するはずです), 衝突処理を行います (弾性衝突またはすり抜けならば速度を交換するだけです).

このあたりから原論文の記述がよくわからなくなるのですが, 次に粒子 $j$ と $j+1$ の衝突時刻 $\tau_j$ をアップデートすれば良さそうです. この結果, おそらく $R$ 分木はヒープ条件を満足しなくなりますから, shiftdown 操作を行ってヒープ条件を満足するようにします. そして, 粒子の情報が更新されたことに伴って $\tau_{j-1}$ と $\tau_{j+1}$ も更新する必要があります. そこで更新後にこれらのデータについて shiftup なり shiftdown なりを行えば再びヒープ条件を満足するようにできます. これらの処理に要する時間は木の深さ程度 $O ( \ln N )$ です.

さて, 再びヒープのルートを見れば次の衝突がいつどこで起こるかがわかります. そこで以上の処理を反復すれば $O ( \ln N )$ で時間発展を追跡することができます.

コメント

原論文のアルゴリズムはほぼこれだけで, あとは具体例に適用した結果などが載っています. 要するに次の衝突時刻をトラックし続けるのにヒープが便利だよ, という話です. 確かにもっともに聞こえますが, いろいろと情報が不足していてこれだけでは十分ではありません.

まず, このアルゴリズムでは情報をアップデートした粒子とアップデートしていない粒子が混在するため, それをどうするか検討する必要があります. 各粒子の情報にそれがいつのものかという時刻を持たせるのが最もわかりやすいでしょうか. 独立時間刻みみたいなものですね.

また, ヒープと粒子情報リストの間で対応関係をすぐに引けるようにしておく必要があります (いちいち探索していたら $O ( N )$ の時間がかかり台無しです). 上の説明ではヒープに粒子 id も保持していますから粒子情報に飛ぶのは簡単ですが, 逆はどうしましょうか. 粒子情報にヒープの中のインデックスも含めておくのが簡単ですかね ($R$ 分木をアップデートする度にそちらも併せて変更せざるを得ないのが面倒ですが).

最後に, ファイル出力する際にすべての粒子の情報をアップデートする必要があり, それには $O ( N )$ の時間がかかるはずです. どの程度の頻度で出力するかによりますが, 交差を詳しく見たいと思っているのならばかなり高頻度で行うことになるでしょうし, 上記アルゴリズムが $O ( \ln N )$ であるとはいえ, そちらの時間が主になりあまりご利益がなさそうです. しかし長時間進化が知りたい場合にはこのアルゴリズムによる時間削減効果はなかなか嬉しいかもしれません.