Position Based Dynamics と Substepping の組み合わせがとても強いという話

最近は The Blob Toy や Pudding など、position based dynamics (PBD) を使った作品をいくつか作っています。

これらの作品ではいずれも通常の陰解法の代わりに substepping を用いているのですが、その効果が想像していたよりもはるかに強力だったので、実際どれほど強いかというのを解説していこうと思います。

Substepping 自体の詳細な説明は元の論文^[1]を読むと分かりやすいです。

陽解法と陰解法

代表的な時間発展の方法の種類に陽解法と陰解法があります。Substepping がどう強力なのか理解するのに必要なので、バネや弾性体の計算について考えたときの両者の特徴について大ざっぱに説明します。

陽解法

陽解法では、あまり難しいことを考えずに素直に速度と位置を公式に従って更新します。例えば計算対象がバネで繋がれた質点なら、現在の自然長との差分から加速度を計算して速度を更新し、その速度を用いて位置を更新する……といった形です。

実際はこの説明だとかなり不正確なのですが、大まかには「現在の時刻の状態から次の時刻の状態を（何らかの式に代入する形で）そのまま計算する」手法だと考えてください。

陽解法のソルバには一般に次のような特徴があります。

• 代入するだけで次の状態が求まるため、高速
• 系のエネルギーが増加する傾向にある
  • 時間刻み幅・バネ定数を大きくするほどこの影響は顕著に表れる
• 不安定になりがち

陽解法の中には高次のルンゲ・クッタ法など精度の高いものも含まれますが、誤差の出方としてはどれも似たような傾向にあります。

陰解法

続いて陰解法です。陰解法では、陽解法と異なり単なる代入によっては次の時刻の値を求めることはできません。代わりに、次の時刻の値が満たすべき方程式を立て、これを解くことで次の時刻の値を求めます。より具体的には、陽解法で現在の時刻の値を用いていた場所を、次の時刻の値に置き換えます^[2]。

すると、次の時刻の値を計算するために次の時刻の値が必要になり、x=f(x) というような形の式が出来上がり、単純な代入では計算できなくなってしまいます。これをどうにかして解くことで次の時刻の値を得るのが陰解法です。

具体的にどうするかは問題設定によって異なり、実は式変形すれば代入で解ける場合や、線形方程式に帰着できる場合、はたまた非線形方程式を解かなければならない場合^[3]などがあります。

陰解法のソルバには一般に次のような特徴があります。

• 方程式を解かなければならないため、低速
• 系のエネルギーが減衰する傾向にある
  • 時間刻み幅・バネ定数を大きくするほどこの影響は顕著に表れる
• 安定する

PBD と陰解法

PBD は基本的には陰解法を用いた手法です。ある時刻の位置から次の時刻における位置に関する（非線形）方程式を立て、それを解いて位置を更新した後に位置の差分から速度を更新する、というのが基本の流れです。

すると、当然 PBD はこの特徴を受け継ぐことになります。

• 系のエネルギーが減衰する傾向にある
  • 時間刻み幅・バネ定数を大きくするほどこの影響は顕著に表れる

すなわち、全体として系のエネルギーは減衰する傾向にあり、しかも時間刻み幅が・バネ定数が大きいほど（すなわちタイプステップあたりの角振動数が大きいほど）エネルギーは速く減衰するのです。

結果として物体を硬くしようとするほどエネルギーが大きく減衰するようになってしまい、振動が ｽﾝｯ…… と素早く収まってしまいます。ひもや布であればこれでも構わないのですが、振動が長く続くような物体を作ろうとすると、ある程度柔らかく設定しなければなりません。

なので、「硬くてプルプルした物体を作るのが難しい」ということになります。

しかし、substepping を使うとこのジレンマを打ち破ることができます。

PBD + substepping

Substepping を用いた PBD では、陰解法で方程式を解くために反復処理を行っている部分を修正し、解を反復修正する代わりにタイムステップを分割します。

つまり、例えば陰解法で 10 回反復していた場合、代わりに一つのタイムステップを 10 分割し、分割されたタイムステップ (substep) では反復を 1 回しか回しません。

反復を 1 回しか回さないということは事実上 substep 内では陽解法が回っているに等しいのですが、PBD に限っては「反復数や時間刻み幅によらず無条件でシミュレーションが安定する」という強い性質があるので、substepping によってエネルギーが発散することはありません。

また、substep 中では高速化のために衝突情報の更新を行いません。相対位置を利用して仮想的な衝突点を補間してやることで、ほとんどの場合問題なく動きます。それでも物体の位置の更新は必要になるので、単純な陰解法に比べると計算量はやや増えることになります。

しかし、それ以上の効果があるのがすごいところです。

Substepping の威力の検証

実際どれほど効果があるのかを簡単なシーンで検証しました。

質点を距離拘束で繋いで作ったひもがあります。

動画は陰解法の反復を 10 回にして計算したものです。理想的には横一直線になるのですが、10 回程度の反復では誤差を収束させ切れず、垂れ下がった状態になります。

ある程度弾性がある（びよんびよんしている）ことにも注目します。

陰解法の反復を増やしたときの様子

反復回数を倍の 20 回にしたものです。先程よりも高い位置で安定していることが分かります。

100 回にするとこうなります。ひもはだいぶ硬くなりましたが、弾性も失われてしまいました。

反復を 1000 回まで増やすと、ひもはほぼ一直線になります。そしてほぼ振動しません。さらに当然ですが、計算時間も 1000 倍かかっています。

以上を踏まえて substepping に切り替えたときの様子を見ていきます。

同じ反復回数で陰解法から substepping に切り替える

トータルの計算回数を 10 回に固定した状態で、完全な陰解法から完全な substepping に段階的に切り替えを行い、挙動がどう変化するかを見ていきます。

比較のために、陰解法で反復を 10 回行った場合の動画を再掲します。

10 反復

次に、2 回の substeps と、各 substep において 5 回の陰解法の反復を行ったものを見てみます。

2 substeps、5 反復

先ほどより安定する位置が少し高くなったのが分かるでしょうか。これはすなわち、拘束がより硬くなったことを意味します。そして、心なしかよく弾むようにもなっているように見えます。

続いて、5 回の substeps で 2 回ずつの陰解法の反復を行ったものを示します。

5 substeps、2 反復

明らかな違いが分かると思います。拘束は硬くなり、よくエネルギーが保存している様子が分かります。

最後に、10 回の substeps で 1 回ずつの陰解法の反復を行ったものがこちらです。

10 substeps

比較のために、100 回の陰解法の反復を行ったものを再掲します。

100 反復

驚くべきことに、たった 10 回の substeps で陰解法の 100 回分に相当する硬さの拘束が得られており、さらに非常によくエネルギーを保存しています。

つまり、substepping を使うと拘束の硬さとエネルギー保存性を両立できるだけでなく、計算の効率も非常に良くなるわけです。同じ硬さの拘束が得たければ、その分ステップの分割数を減らして計算を高速化させることができます。

ちなみに、10 回の substeps で 100 回分の反復に相当する硬さになったのは偶然ではなく、一般に N 回の substeps で概ね N^2 回分の反復の硬さの拘束を得ることができます。これはなかなかに驚異的です。

Substeps をさらに増やすとどうなるか

ステップの分割数をさらに増やしたときの挙動を見ていきます。

20 substeps では次のようになります。

100 substeps だとこうなります。

これは陰解法ではおおよそ 10000 回の反復に相当する硬さになっており、ひもは一切たわんでいないように見えます。

そこで、初期変位を加えてシミュレーションを実行してみると次のようになります。

ギターの弦のように振動している様子が確認できます。ステップを分割しているため、本来の時間刻み幅では捉え切れなかった挙動まできちんと再現することができます。

例えば、60 FPS のシミュレーションではナイキスト周波数の関係で 30 Hz の振動までしか原理的に捉えることができませんが、100 回もの substeps があれば最大 3000 Hz の振動まで捉えることができます。

なぜここまで硬い拘束を実現できるのか？

こうなるとなぜここまで substepping が優秀なのかというのが疑問になってきます。自分なりに考えてみた結果、「勢いが乗るから」ではないかという結論が得られました。

陰解法の反復というのは、基本的に現在の暫定的な解を一定の割合で真の解に近付ける効果を持ちます。例えば、一度の反復で真の解に 30% 近付くと仮定した場合、反復を 3 回繰り返すと誤差は 70% → 49% → 34% という風に減っていきます（1 次収束といいます）。

このとき、それぞれの誤差の減少は各反復で独立しています。どういうことかというと、1 回目の反復で減った誤差というのは、2 回目の反復においては特に役立っていないということです。まあ当然の話で、N 回目の反復で減る誤差というのは N 回目の反復で行った計算によるもののみであるということです。N – 1 回目以前の反復で何が起きたかというのは、初期状態の変化を除いて N 回目の反復において影響を及ぼしません。

しかし、substeps を考える場合は少し話が変わってきます。各反復で得られた修正量が速度に乗るのです。1 substep 目で修正された位置は（位置ベース物理の計算式に従い）速度に上乗せされ、慣性の法則によって 2 substep 目に引き継がれます。2 substep 目での修正量は同じく速度に乗るので、3 substep 目では 1 と 2 双方の substep での影響を受けることになります。同様に、4 substep 目では 1, 2, 3 全ての substep での影響を受けることになります。

結果として、N substep 目では 1, 2, …, N – 1 substep 目全ての影響を受けることになり、N 回の substep での計算を終えた時点で通常の反復の 2 乗に相当するオーダーの修正^[4]が実現できているのではないかと思います。さらに、この「勢い」というのはエネルギーそのものなので、エネルギーの保存性と相乗効果があるということになります。

結果的に、The Blob Toy のシミュレーションは 10 回の substepping でこれだけの安定性とエネルギー保存性を確保できています。

ちなみに通常の PBD としての反復 10 回だとこうなります。圧倒的な違いがあることが分かります。

Pudding では 8 回の substepping でこれだけ安定したシミュレーションが実現できています。これも通常の 8 反復ではありえない話です。

位置ベースでないシミュレーションへの substepping の利用

これまで説明した通り substepping は特に PBD において高い効果を発揮するのですが、位置ベースでないシミュレーションにおいても、細かい substep への分割は計算の非線形性を大きく緩和する効果があるので、うまく利用すれば計算効率を上げることが期待できます。

物理エンジンは一般に非線形性が苦手、というのはだいぶ前に書いた記事の通りです。

とはいえ現在普及している state of the art を substepping で更新するのはなかなか難しい……と思っていたのですが、なんとあの剛体物理エンジン Box2D がやってくれました。

つい最近リリースされた Box2D 3.0 では、substepping を利用したソルバがデフォルトで搭載されています。というか多分他のソルバのオプションがなさそうで、これは全面的に新しいソルバが既存のソルバの性能を上回ったということを意味します^[5]。

実際は substepping だけでなく soft constraint のテクニックも併用されており、これらの併用によって剛体計算における substepping の効果を最大限に発揮できたようです。

Box2D 3.0 で利用されているソルバの詳細については、作者の Erin Catto 氏による Solver2D についての記事に書かれています。とても興味深い内容なので気になる方は是非動画も併せて見てみてください。この辺りの開発力というか研究力は本当にさすがだと思います。

いろんなソルバによる様々なシーンでの性能を比較した動画です。新ソルバの前身である TGS_Soft の威力が分かります。