PBDとは｜Position Based Dynamicsの仕組みとXPBDの違いを実装コードで解説

PBD（Position Based Dynamics）は、力から加速度を求めて積分する代わりに、質点の位置を直接補正して拘束条件を満たす物理シミュレーション手法です。Matthias Müllerらが2006年のワークショップVRIPHYSで発表し、布・ロープ・髪・ソフトボディをゲームで安定して動かす方法として広まりました。後継のXPBD（Extended Position Based Dynamics）は、PBDの弱点だった「硬さが反復回数と時間刻みで変わる」問題を解いた拡張です。

この記事では、PBDの1ステップの処理手順、距離拘束の式、XPBDとの違いをMüllerらのPBD原論文とXPBD論文に沿って説明し、Node.jsで実行を確認した約40行のコードと、その実測値で両者の差を示します。なお「PbD」はプライバシー・バイ・デザイン、「.pbd」は別形式のファイル拡張子、「PBD」は企業の略称としても使われますが、この記事が扱うのは物理シミュレーションのPBDだけです。

まとめ：PBDとXPBDの要点

• PBDの考え方：速度で位置を予測し、拘束（距離・曲げ・衝突など）を満たすように予測位置を反復で補正し、補正後の位置から速度を逆算する。
• 強み：時間刻みを大きくしても発散しにくく、衝突は「めり込んだ点を外へ押し出す」だけで扱える。実装が短い。
• 弱点：拘束の硬さ k が反復回数と時間刻みに依存する。後述の実測では、同じ k=0.5 の鎖が反復5回で30.56%、80回で1.76%伸びた。
• XPBD：拘束ごとにコンプライアンス α（剛性の逆数）とラグランジュ乗数 λ を持たせ、材料の硬さを α で指定できるようにした。後述の鎖では、反復20回と80回の静止時の伸びが19.64%と19.62%でほぼ一致した（反復が少なく収束しきらない場合の誤差は残る）。
• 今から実装するなら：PBDでなくXPBDを選び、反復回数を増やすよりサブステップ（1フレームを細かく割る）を増やす。
• 本記事のサンプルの限界：距離拘束だけの鎖は挙動を確かめるためのもので、応力やひずみの定量評価には使えない。定量評価には材料モデルとメッシュを持つ有限要素法の定式化と、実験値との照合が要る。

Position Based Dynamicsの定義と力ベース手法との違い

一般的な物理シミュレーションは力ベースです。ばねの力や重力を合計し、ニュートンの運動方程式で加速度を求め、速度、位置の順に積分します。この方式で硬いばねを陽解法で解くと、時間刻みが大きいときに補正が行き過ぎて振動が増え、最後は発散します。

PBDは速度の層を飛ばして位置を直接操作します。原論文の要旨は、剛体シミュレータの多くが速度を直接操作する（インパルスベース）のに対し、PBDは「速度の層も省いて位置に直接作用する」と説明し、主な利点に制御しやすさと、陽的積分の行き過ぎ（overshooting）を避けられる点を挙げています。著者4人は当時PhysXを開発していたAGEIAの所属で、論文ではこの手法でゲーム向け物理ライブラリのリアルタイム布シミュレータを作ったと述べています。

論文は積分方式の位置づけも書いています。反復1回なら陽解法に近く、反復を増やすと拘束系をいくらでも硬くでき、陰解法に近い振る舞いになります。ボトルネックは衝突判定からソルバーへ移ります。

PBDのアルゴリズム：1ステップの処理手順

原論文のアルゴリズムを、1フレーム分の処理として並べると次のとおりです。w は質量の逆数で、w=0 の点は固定点になります。

1. 外力（重力など）で速度を更新する：v ← v + Δt·w·f_ext
2. 速度を減衰させる（dampVelocities）
3. 予測位置を求める：p ← x + Δt·v（陽的オイラーの1ステップ）
4. 現在位置 x から予測位置 p への移動を調べ、衝突拘束を生成する
5. すべての拘束について、予測位置 p を補正する処理を決められた回数だけ繰り返す（solverIterations）
6. 補正後の位置から速度を逆算する：v ← (p − x) / Δt
7. 位置を確定する：x ← p
8. 摩擦や反発を速度に反映する（velocityUpdate）

要は手順3・5・6・7の4つです。論文によると、6と7は陽解法のように未来へ外挿するのではなく、ソルバーが求めた「物理的に妥当な位置」へ点を移すため、不安定になりうるのはソルバー自身だけで、その安定性は時間刻みではなく拘束関数の形で決まります。

衝突の扱いもこの流れに乗ります。手順4で「床より下に入ったら床の高さへ戻す」といった不等式拘束を毎ステップ作り直し、手順5で他の拘束と一緒に解きます。めり込んだ点を正しい位置へ射影するだけで貫通を完全に解消できる点は、論文の要旨でもPBDの利点として挙げられています。

距離拘束の定式化と制約投影

布やロープの基本になる距離拘束を例に、手順5の中身を見ます。2点 p₁、p₂ の距離を静止長 d に保つ拘束は、次の式で表します。

C(p₁, p₂) = |p₁ − p₂| − d = 0

C が0でなければ、2点を結ぶ方向 n = (p₁ − p₂) / |p₁ − p₂| に沿って、質量の逆数の比で補正量を配分します。

• Δp₁ = −w₁ / (w₁ + w₂) · C · n
• Δp₂ = +w₂ / (w₁ + w₂) · C · n

重い点（w が小さい点）ほど動かず、固定点（w=0）は動きません。どちらの点も有限の質量を持つ場合、質量で重み付けした補正量の和 m₁Δp₁ + m₂Δp₂ は0になるので、この内部拘束だけで2点の重心が動くことはありません（原論文も運動量が保存されることを示しています）。固定点や外部との衝突拘束が絡む場合はこの限りではありません。曲げ・体積・衝突の各拘束も、拘束関数 C とその勾配を定義すれば同じ枠組みで扱えます。

ガウス・ザイデル型とヤコビ型の使い分け

原論文のソルバーは、拘束を1本ずつ解いてすぐ位置に反映するガウス・ザイデル型です。前の拘束の補正が次の拘束に即座に効くので収束は速い一方、拘束の処理順に依存するため並列化しにくくなります。GPUで大量の拘束を並列に解く実装では、全拘束の補正量をいったん粒子ごとに集め、平均してから反映するヤコビ型が使われます。後述のFleXの元になった2014年の論文がこの方式で、ヤコビ型は条件によって収束が保証されないため、平均化で補正量を抑えて安定させています。GPUでもガウス・ザイデル型を使いたい場合は、点を共有しない拘束どうしをグループ（色）に分け、グループ内だけを並列に解く方法があります。収束の速さを優先するならガウス・ザイデル型、実装の単純さと並列度を優先するならヤコビ型、という軸で選びます。

実装例：JavaScriptで鎖をPBDとXPBDで解く

次のコードは、長さ5cmの区間19本でつないだ質点20個の鎖を上端で固定して吊るし、600フレーム（60fpsで10秒）後の全長が静止長から何%伸びたかを出力します。mode=’pbd’ は硬さ k で補正量を縮める原論文の方式、mode=’xpbd’ はコンプライアンス α を使うXPBDの方式です。Node.js v26.5.0でそのまま実行できます。

// 鎖（質点N個）を上端固定で吊るし、静止後の全長の伸びを測る
function simulate({ n = 20, iters = 10, substeps = 1, mode = 'pbd', k = 0.5, alpha = 1e-4, steps = 600 }) {
  const dt = 1 / 60, h = dt / substeps, L = 0.05, g = -9.8, damp = Math.pow(0.99, 1 / substeps);
  const x = [], y = [], px = [], py = [], vx = [], vy = [], w = [], lambda = [];
  for (let i = 0; i < n; i++) { x.push(i * L); y.push(0); px.push(0); py.push(0); vx.push(0); vy.push(0); w.push(i === 0 ? 0 : 1); }
  for (let s = 0; s < steps * substeps; s++) {
    for (let i = 0; i < n; i++) { if (w[i] > 0) vy[i] += g * h; px[i] = x[i] + vx[i] * h; py[i] = y[i] + vy[i] * h; }
    for (let j = 0; j < n - 1; j++) lambda[j] = 0;
    for (let it = 0; it < iters; it++) {
      for (let j = 0; j < n - 1; j++) {
        const a = j, b = j + 1, wSum = w[a] + w[b];
        if (wSum === 0) continue;
        const dx = px[a] - px[b], dy = py[a] - py[b], d = Math.hypot(dx, dy);
        if (d === 0) throw new Error(`particles ${a} and ${b} overlap`);
        const C = d - L, nx = dx / d, ny = dy / d;
        let dl;
        if (mode === 'pbd') dl = -k * C / wSum;
        else { const at = alpha / (h * h); dl = (-C - at * lambda[j]) / (wSum + at); lambda[j] += dl; }
        px[a] += w[a] * dl * nx; py[a] += w[a] * dl * ny;
        px[b] -= w[b] * dl * nx; py[b] -= w[b] * dl * ny;
      }
    }
    for (let i = 0; i < n; i++) { vx[i] = (px[i] - x[i]) / h * damp; vy[i] = (py[i] - y[i]) / h * damp; x[i] = px[i]; y[i] = py[i]; }
  }
  let len = 0;
  for (let j = 0; j < n - 1; j++) len += Math.hypot(x[j] - x[j + 1], y[j] - y[j + 1]);
  return ((len / ((n - 1) * L) - 1) * 100).toFixed(2) + '%';
}
for (const iters of [5, 20, 80]) {
  console.log(`iters=${iters}  PBD(k=0.5)=${simulate({ iters, mode: 'pbd' })}  XPBD(alpha=1e-4)=${simulate({ iters, mode: 'xpbd' })}`);
}
console.log(`substeps=20,iters=1  XPBD=${simulate({ iters: 1, substeps: 20, mode: 'xpbd' })}  PBD=${simulate({ iters: 1, substeps: 20, mode: 'pbd' })}`);
console.log(`substeps=1,iters=20  XPBD=${simulate({ iters: 20, substeps: 1, mode: 'xpbd' })}`);

PBDとXPBDの差は、内側ループの dl を求める1行だけです。質点の質量は1kg（w=1）、長さの単位はm、コンプライアンス α の単位はm/N です。速度に掛けている damp は振動を止めるための簡易な減衰で、原論文の手順2（dampVelocities）の代わりです。サブステップ数を変えても1フレームあたりの減衰が0.99倍でそろうよう、0.99 のサブステップ数乗根にしています。2点が完全に重なると方向を定義できないため、その場合は例外を投げます。実行結果は次のとおりでした。

iters=5  PBD(k=0.5)=30.56%  XPBD(alpha=1e-4)=24.04%
iters=20  PBD(k=0.5)=7.49%  XPBD(alpha=1e-4)=19.64%
iters=80  PBD(k=0.5)=1.76%  XPBD(alpha=1e-4)=19.62%
substeps=20,iters=1  XPBD=19.62%  PBD=0.38%
substeps=1,iters=20  XPBD=19.64%

PBDの弱点：硬さが反復回数と時間刻みで変わる問題

上の結果のPBD列が、この手法の最大の弱点です。k=0.5 という同じ設定の鎖が、反復5回では30.56%伸び、20回で7.49%、80回で1.76%まで硬くなりました。1フレームを20サブステップに割って反復1回ずつにすると0.38%で、ほぼ伸びないロープになります。反復回数や時間刻みを変えるたびに、アーティストが決めた「布の硬さ」が別物になるということです。

原論文もこの問題を認めています。k を補正量に掛けるだけだと、1本の距離拘束の残差は ns 回の反復後に Δp(1 − k)^ns となり、効き方が反復回数に対して非線形になります。そこで論文は k の代わりに k′ = 1 − (1 − k)^1/ns を掛け、残差を Δp(1 − k) に揃える補正を示しました。ただし同じ段落で、この補正後も材料の硬さは時間刻みに依存したままであり、固定時間刻みのリアルタイム環境なら問題にならない、と限定しています。

ゲームのように時間刻みも反復回数も固定で、見た目が破綻しなければよい用途なら、PBDのままでも困りません。硬さをパラメータとして資産に保存し、品質設定で反復回数を変えるような作り方をするなら、次のXPBDを使うべきです。

XPBDとは：コンプライアンスで硬さを反復回数から切り離す拡張

XPBDは、Miles Macklin、Matthias Müller、Nuttapong Chentanez（いずれもNVIDIA）が2016年のMotion in Games（MIG ’16）で発表した「XPBD: Position-Based Simulation of Compliant Constrained Dynamics」で提案されました。論文は要旨で、PBDの長年の問題である「反復回数と時間刻みに依存する拘束の硬さ」に取り組むと書いています。

論文のAlgorithm 1は、元のPBDに3行（λ の初期化、Δλ の計算、λ の更新）を足しただけの形をしています。拘束ごとにコンプライアンス α（剛性の逆数）を持たせて α̃ = α / Δt² とし、拘束ごとの累積ラグランジュ乗数 λ を各ステップの初めに0へ戻してから積み上げます。1本の拘束を解くたびに求める増分は次の式です。

Δλ_j = (−C_j − α̃_jλ_j) / (∇C_j M⁻¹ ∇C_j^T + α̃_j)

位置の補正は Δx = M⁻¹∇C^TΔλ です。論文が指摘するとおり、α=0（無限に硬い拘束）のときこの式は、硬さ係数 k=1 のPBDの補正と一致します。既存のPBD実装に数行足すだけで移行できるのはこのためです。上のコードの else 節がこの式そのもので、距離拘束では ∇C M⁻¹ ∇C^T が w₁ + w₂ になります。

実測では、α=1e-4 の鎖の伸びは反復20回で19.64%、80回で19.62%、20サブステップ×反復1回で19.62%でした。区間 j にかかる張力はその下にぶら下がる質点の重さなので、α×張力の総和から静止時の伸びを手計算すると 1e-4 × 9.8 × (1+2+…+19) = 0.186m、静止長0.95mに対して19.6%です。反復20回以上とサブステップ20回の結果はこの理論値とよく一致しており、少なくともこの鎖の静止時の伸びについては、反復回数や時間刻みを変えても α が決める値に落ち着くことが確認できます。反復5回の24.04%だけが外れているのは、反復が少なすぎて収束しきっていないためです。α は材料の硬さを反復回数から切り離しますが、収束しきらない分の数値誤差まで消すわけではありません。

XPBDの計算量の配分：反復回数よりサブステップ数

XPBDの計算量の配分については、2019年のSCAで発表された「Small Steps in Physics Simulation」（Macklinほか、NVIDIA）が指針を出しています。要旨は、1回の大きな時間刻みでソルバーを n 回反復するより、時間刻みを n 分割して各サブステップで反復を1回だけ行うほうが効果が高い、というものです。論文はこの方式を、陰解法の大きな時間刻みと比べて拘束誤差と数値減衰が大幅に小さく、陽解法より広い剛性範囲で安定すると報告しています。

サブステップごとに衝突判定をやり直すと費用が膨らむため、同論文のAlgorithm 1は衝突判定をフレームの先頭で1回だけ行い、サブステップ内ではその結果の拘束を解く構成です。論文は限界も書いています。Δt² の項があるため、サブステップを重ねると32ビット浮動小数点の精度の限界に達することがあり、大きな座標値に小さな位置の変化を足しても値が変わらない場合がある、というものです。論文の例はすべて32ビットで、反復をさらに増やすなら倍精度が必要になりうるとしています。原点から遠い座標で細かいステップを刻むときは、浮動小数点数の仕組みと誤差を踏まえて座標系の置き方を決めてください。

PBDの応用：布・髪・ソフトボディ・流体・剛体

布・ロープ・髪の距離拘束と曲げ拘束

距離拘束で布の伸びを、隣り合う三角形の二面角を保つ曲げ拘束で折れにくさを表します。Müllerらの2006年の原論文も、この2種類の拘束と自己衝突で布を実装しています。髪やロープは質点の列に距離拘束と曲げ拘束を並べたもので、上の鎖のコードがそのまま出発点になります。

流体の密度拘束とPosition Based Fluids

MacklinとMüller（NVIDIA）がSIGGRAPH 2013で発表し、ACM Transactions on Graphics 32(4)に掲載された「Position Based Fluids」（PBF）は、PBDの枠組みに密度の拘束を加えた流体の手法です。粒子 i ごとに C_i = ρ_i / ρ₀ − 1（ρ₀ は静止密度、ρ_i は近傍粒子からSPHのカーネルで推定した密度）という拘束を置き、密度を一定に保つことで非圧縮性を表します。論文の要旨は、現代のSPHソルバーに近い非圧縮性と収束性を持ちながら、PBDの安定性を引き継いで大きな時間刻みを使え、リアルタイム用途に向くと述べています。

剛体・関節のXPBDによる処理

剛体は2020年のSCAで発表された「Detailed Rigid Body Simulation with Extended Position Based Dynamics」（Müllerほか、NVIDIA）でXPBDに取り込まれました。論文は、接触と複数種類の関節、ソフトボディとの相互作用を扱う剛体ソルバーを実装するための詳細をすべて示したとしています。速度レベルで拘束を線形化する従来の剛体シミュレータと違い、常に最新の拘束方向で解くため、曲面に高速でぶつかる物体も追跡しやすくなります。前述のSmall Stepsのサブステップ方式と組み合わせる前提の設計です。

PBD・XPBDを使える主なツールとライブラリ

映像制作ではHoudiniのVellumが代表例です。SideFXの公式ドキュメントは、Vellumを「拡張Position Based Dynamics（XPBD）のアプローチを使うシミュレーションの枠組み」と説明し、布・髪・ソフトボディ・風船・粒を作れるとしています。Houdiniの特徴とBlenderとの違いを把握したうえで、布や髪の表現にVellumを選ぶのが近道です。

リアルタイム用途では、PhysX 5に PxPBDParticleSystem というクラスがあります。PhysX 5.1.3の公式ドキュメントは、この粒子システムのソルバーが粒子間の力学をPBDで解き、流体や変形体を扱えると説明しています。ただし布の扱いは版によって変わっています。NVIDIAのOmniverse Physicsの移行ガイドは、旧実装の粒子布（particle cloth）がPBDの質点ばねモデルだったのに対し、移行先の面の変形体（surface deformables）はXPBDによる有限要素法（共回転線形弾性）で、物理モデルが違うため同じ挙動は近似しかできないと書いています。PhysXを物理演算に使うIsaac Simのような環境で布を扱うときは、使っている版がどちらのモデルかを先に確認してください。XPBDと有限要素法が組み合わせて使われている例でもあります。

FleXは「Unified Particle Physics for Real-Time Applications」（SIGGRAPH 2014）の系譜にある、すべての物体を粒子で表すリアルタイム物理ライブラリです。この論文は、GPUで並列に解くために拘束をヤコビ型で解き、粒子ごとに補正量を平均する方式（constraint averaging）をPBDに導入しています。GitHubのリポジトリは2021年4月から更新が止まっているため、新規採用より手法の参考資料として読むのが妥当です。自分で実装して中身を理解したいなら、Jan BenderらのC++ライブラリ「PositionBasedDynamics」（MITライセンス）が剛体・変形体・流体の実装例を一通り含んでいます。最新リリースは2022年12月の2.2.0ですが、リポジトリへのコミットは2026年9月時点も続いています。

有限要素法・SPH・質点ばねモデルとの比較と採用しない場面

表の5つは同じ層の選択肢ではありません。有限要素法やSPHは物体をどう離散化するかの方法で、PBDやXPBDは拘束をどう解くかの方法です。そのため前述のOmniverse PhysicsのようにXPBDで有限要素法の弾性モデルを解く組み合わせや、SPHの密度推定を使うPBFのような組み合わせが成り立ちます。

質点ばねモデルは、ばね定数を上げるほど安定に必要な時間刻みが小さくなります。PBDはこの制約から自由になった代わりに、硬さの物理的な意味を失いました。XPBDはその意味をコンプライアンスとして取り戻しましたが、論文自身が「陰的な運動方程式を厳密に解いているとは言えない」と書いているとおり、近似解法です。

判断を分けるのは、結果の数値を保証する必要があるかどうかです。部材の応力や変形量を設計値として使う検証に、本記事のような距離拘束だけのモデルを使うべきではありません。材料モデルとメッシュ、収束判定、実験値との照合までそろえたCAEで使われる有限要素法の解析が必要です。XPBDを使う場合でも、拘束を有限要素法の弾性エネルギーから作り、収束するまで反復して誤差を評価しない限り、設計値の根拠にはなりません。反対に、60fpsなら1フレーム約16.7ms以内に結果を返す必要があり、多少の誤差より破綻しないことが優先されるゲーム・VR・インタラクティブな演出では、XPBDが第一候補になります。流体では、PBF自体がSPHの密度推定を使い、原論文も既存のSPHソルバーと非圧縮性・収束性を比べています。PBFを選ぶ理由は、大きな時間刻みで安定させたい場合と、布や剛体と同じ拘束ソルバーで流体を連成させたい場合です。

よくある質問

PBDとXPBDの違いは何ですか？

XPBDは拘束ごとにコンプライアンス α と累積ラグランジュ乗数 λ を追加したPBDの拡張です。PBDでは設定した硬さの効き方が反復回数と時間刻みで変わりますが、XPBDでは材料の硬さを α で指定でき、反復が収束していれば結果は α で決まる値に近づきます。α=0 のXPBDは、硬さ係数 k=1 のPBDと同じ補正になります。

PBDは何の略ですか？

物理シミュレーションの文脈ではPosition Based Dynamics（位置ベース動力学）の略です。同じ綴りでも、個人情報保護のPrivacy by Design（PbD）や企業名の略称、「.pbd」という別形式のファイル拡張子を指す場合があり、これらはPosition Based Dynamicsとは無関係です。

PBFとは何ですか？

Position Based Fluidsの略で、PBDの枠組みに粒子の密度を一定に保つ拘束を加えて非圧縮性の流体を表す手法です。MacklinとMüllerが2013年に発表しました。

反復回数とサブステップ数はどちらを増やすべきですか？

XPBDではサブステップ数を増やすほうが効果的です。「Small Steps in Physics Simulation」（2019年）は、1回の時間刻みでn回反復するよりn個のサブステップで1回ずつ反復するほうが拘束誤差と数値減衰が小さいと報告しています。

PBDのソルバーはGPUで並列化できますか？

できます。代表的な方法は2つあり、各拘束の補正量をまとめて計算し粒子ごとに平均して反映するヤコビ型と、点を共有しない拘束をグループに分けてグループ内を並列に解くガウス・ザイデル型です。ヤコビ型は収束が遅くなりやすいので、反復回数かサブステップ数で補います。

関連記事

• CAEとは？読み方と定義・解析種別の分類とCAD/CAMの役割分担を解説
• Houdiniとは？何ができるソフトか・特徴とBlenderとの違いをわかりやすく解説
• Isaac Simの基本概要とNVIDIAロボティクス開発で果たす役割
• 浮動小数点数とは？IEEE 754の仕組みと誤差・型の選び方を実装目線で解説