計算機圖形中動態系統模擬最流行的方法是基於力的。累積內部和外部力量，根據牛頓的第二個運動定律計算加速度。然後使用時間積分方法來更新速度，最後是物件的位置。
一些模擬方法（大多數剛性體模擬器）使用基於衝量的方法並直接操縱速度。
PBD是一種省略了速度層的直接作用於位置的控制方法，方法的主要優點是它的可控性、效能、數值穩定，在遊戲開發領域應用很廣，近年的UE chaos底層也採用了PBD方法。

Verlet積分

verlet積分在遊戲開發領域應用也很廣泛，並且PBD論文中提到了PBD的思想和verlet有相近之處，這裡也順便介紹一下。

基本積分方法

積分方法的階數一般是看泰勒展開中導數的最高項。n階方法的誤差是\(O(h^{n+1})\)。
值得一提的是，explicit euler、semi-implicit euler都是一階方法。
積分方法討論的通常是給定一個函數S(t)，然後我們沒有辦法得到S(t)的解析式（假定這裡是和時間有關的函數），此時我們通過一些方法來對S(t) 做近似，達到一個我們給定一個時刻的資料，可以推算其他時刻資料的效果，這裡一切的理論基礎都是泰勒展開，然後進行各種化簡。
這裡用圖來表示，小紅色的矩形代表積分過程，左邊是explicit euler，右邊是semi-implicit euler。

左圖用的矩形的左邊來表示一個小時間步裡整體的情況，意思就是時間步剛開始的值。左圖用的矩形的右邊來表示一個小時間步裡整體的情況，意思就是時間步結束時的值。

Verlet 算位置

verlet積分是一種牛頓運動方程的積分方法，他的積分方法有著良好的數值穩定性，對比簡單的尤拉方法沒有顯著的額外計算成本。

\[ x _{t+\Delta t} = x_t + v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 + \frac{1}{6}b(t)\Delta t^3 + O(\Delta t^4) \]

\[ x _{t-\Delta t} = x_t - v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 - \frac{1}{6}b(t)\Delta t^3 + O(\Delta t^4) \]

兩者求和，消去一階和三階項得到

\[ x _{t+\Delta t} = 2 x_t - x _{t-\Delta t} + a(t)\Delta t^2 + O(\Delta t^4) \]

我們可以將\(t+\Delta t\) 當成當前幀所在的時刻，那麼t就是上一幀，\(t-\Delta t\)就是上上一幀，我們只需要兩個位置資訊和上一幀的力資訊，就可以得到當前幀的位置資訊，並且這個位置資訊是三階精度的，誤差四階精度，wiki百科中提到，這個方法這裡的誤差是區域性誤差，並且公式中的加速度是根據精確解計算的，全域性誤差（即累積誤差）與時間步長的平方（Δt^2）成正比。雖然區域性離散化誤差為 4 階，但由於微分方程的2階，推導一番後全域性誤差為 2 階，我們將它定義為2階方法。

有的書上也說這個公式也可以從二階中心差分推匯出，所以是二階方法。

Verlet 算速度

接下來是求速度資訊，
兩者求差，消去二階項得到

\[ x(t+\Delta t)−x(t−\Delta t)=2v(t)\Delta t+ \frac{1}{3}b(t)\Delta t^3 \]

這樣可以推算上一幀的速度資訊

\[ v(t) = \frac {x(t+ \Delta t) - x(t - \Delta t)}{2\Delta t} + \frac{b(t)\Delta t^3}{3\Delta t} \]

\[v(t) = \frac {x(t+ \Delta t) - x(t - \Delta t)}{2\Delta t} + O(\Delta t^2) \]

這個算上一幀速度的方法是一階精度二階誤差的。
在他的原始方法裡推算當幀的速度用的中值定理+近似，導致速度一定是不準的，得出速度公式為

\[v(t+\Delta t) = \frac {x(t+ \Delta t) - x(t)}{\Delta t} + O(\Delta t) \]

這個速度的準度可以說幾乎沒有，零階精度，一階誤差，推導可以看verlet一開始 \(x(t+\Delta t)\)的展開將xt移到左邊再除以步長得到：

\[ \frac {x(t+\Delta t)−x(t)}{\Delta t} =v(t)+ \frac{1}{2}a(t)\Delta t + \frac{1}{6}b(t)\Delta t^2 \]

所以如果\(\frac{1}{2}a(t)\Delta t + \frac{1}{6}b(t)\Delta t^2\)作為誤差，誤差就是\(O(\Delta t)\)的。

wiki百科裡則提供了一種更加高精度的速度推算方法叫 velocity verlet ，利用梯形法（頭尾值加起來求平均，之前的小矩形可以想象成梯形來看）算的加速度來推算速度，讓速度的精度到達二階，至於為什麼是二階，因為梯形法的結果是二階精度，這個精度貌似可以傳遞。

\[v(t+\Delta t) = v(t) + \frac {a(t)+a(t+ \Delta t)}{2} \Delta t \]

a的\(t+\Delta t\)是一般來說是一個和位置有關的函數，先求了位置，接下來這裡面都是已知量了。

PBD

基於力的方法解碰撞

首先一個直觀印象關於如何處理動力學中的穿插問題，基於力的方法往往是計算穿插深度，然後給根據深度和穿插的方向給予力，這過程中有一個剛度係數，表示這個物體的堅硬程度，調的越大，兩個物體越難穿插，你就會看到物體看起來很硬，實際航就是穿插產生的力的縮放係數，剛度越大實際上要求時間步長越小，否則很容易出現穿插過深彈出一個很大的力這種情況，其實就是數值解誤差太大。

過沖問題

顯式積分法，也稱為開環法或單步法，廣泛用於求解科學和工程各個領域的常微分方程 (ODE)。然而，這些方法可能會遇到過沖問題，這會導致不準確或不穩定的解決方案。
產生原因：
1、兩剛體相交過深導致計算出來的力太大，求解的初始條件殘差太大
2、時間步長太大
3、高頻振盪：在求解涉及高頻振盪的 ODE 時，顯式方法可能難以準確捕獲這些振盪，從而導致超調。在這種情況下，使用對剛性系統更穩定的隱式積分方法會有所幫助。
4、誤差累積：使用顯式方法時，誤差會隨著時間的推移而累積，從而導致顯著的超調。

要解決顯式積分中的過沖問題，請考慮使用更小的時間步長、自適應時間步長、高階方法，或者在適當的時候切換到隱式或半隱式方法。此外，優化初始條件並使用靈敏度分析可以提高解決方案的準確性和穩定性。

基於位置的方法解碰撞

在PBD方法中，當檢測到兩個物體發生穿透時，直接根據約束脩正物體位置到一個不會碰撞的位置，然後更新速度資訊。他其實是一個反向的過程，雖然這中間力的計算不明確了，但是表現是正常的，這是我們期待的，這就是pbd論文提到的visually plausible 視覺正確，並且避免了之前提到的過沖問題。

演演算法流程

我們表達一個物體，使用N個頂點，M個約束，一個頂點i對應的質量是mi，位置則是xi，速度 vi
約束的索引一般使用j
我們對於其中一個約束有：

定義一個約束影響的頂點數量是nj，可以理解為第 j個約束所影響的頂點數目為nj
定義約束函數C 頂點數目為3nj，位置資訊為xyz，所有的資料總量就是3nj，約束函數就是讀取這個3nj的資料算出一個實數
約束函數影響的頂點的索引，nj個索引
定義這個約束的剛度值
定義這個約束是等式約束和不等式約束
等式約束比較強硬，要滿足約束函數等於0，不等式約束則是約束函數大於等於0就行，kj定義了這個約束的剛度，剛度越大這個約束要越快被滿足。

演演算法最開始是先初始化x和v，還有質量的倒數
每個迴圈開始，先通過力算速度，在通過算出的速度算位置，計算出一個預測的速度、位置。
然後在根據這個預測的位置情況進行碰撞約束的生成。
然後走一個固定的迭代次數，做約束投影，讓p，其實就是位置，落在正確的滿足約束函數的地方，因為有多個約束，其實這裡就是一個多個約束方程的方程組。
約束投影完了只有得到正確的位置，然後根據投影后的位置和上一次迭代的位置算一下速度（這裡和verlet簡單的算速度方法很像），然後把投影后的位置更新到位置資料裡。
最後在(16)行根據摩擦(friction)和恢復(restitution)係數修改速度。

這個方案是無條件穩定的，因為他不會像一般的explicit方法那樣對未來進行外插，而是將頂點移動到正確的位置。
他的穩定性不取決於時間步長，而是取決於約束函數的樣子。
這個求解器並不像以往的顯式或者隱式積分器，而是處於一種中間態，如果每個時間步只跑一個迭代他看起來像explict，而新增迭代次數，可以更像隱式積分方案。

求解器借用的思想

需要注意的是約束函數是通常都是非線性的，比如距離約束

\[C = |p_1-p_2| - d \]

這裡的非線性，其實指的是他這個方程組的未知數不是一次的，有絕對值的話，可以理解成平方再開方，一般的迭代法解線性方程組使用jacobi、gauss-seidel迭代方法，他們的未知數有多個，但是都是一次的。
但是在這裡作者是借用了GS方法的思想，來解一個非線性方程組。
原gauss-seidel迭代法是，我們一整個方程組，然後從第一條方程開始，代入一個初值，然後方程左邊算出第一個未知數X，然後這個X算出來後會直接替換掉下面方程組求解時用到X的地方，第二條方程算出Y，也是一樣的，先滿足一條方程，然後將它的影響帶到下一條方程，經過固定次數的迭代後，系統的未知數們會收斂，就求解成功了。

值得注意的是，GS方案比較依賴方程的組的滿足順序，過度約束且沒有保證求解順序會導致結果振盪。

關於動量守恆

PBD方法裡提到了動量守恆的話題，PBD本身方法的流程他是想說明按照他的流程走，只要約束函數定義時候是動量守恆的，那麼結果是動量守恆的。
如果你這個系統最後輸出的結果是沒有違背動量守恆的基礎，那麼將看到一些不期待的運動，這裡坐著定義為 Ghost force，看著就像有力量拖著物體旋轉物體一樣。

後面作者提到，對於一個系統來說只有內部約束需要考慮動量守恆，因為外力作用會作用在全域性而產生動量。內部約束不會改變剛體運動狀態，即剛體的旋轉、位移，作者稱之為與剛體模態正交。
因此我們沿著約束函數梯度的方向去影響整個系統他最後滿足約束了約束函數（滿足這個內部力的性質），他的最後就是動量守恆的。
這裡我需要說的詳細的一點，因為我自己也對動量守恆這塊比較困惑：
我的理解是PBD本身其實和動量守恆不沾邊，主要是約束函數的定義是動量守恆的，你沿著約束函數的梯度去影響系統最後達到了約束函數，那麼約束函數的最終狀態是動量守恆的，我們拿這個狀態去做表現，表現就是動量守恆的，但是如果我們因為很多約束，最後求解沒成功，有很多約束沒滿足，那麼最後的表現就是不動量守恆的。
如果我們定義了一個動量不守恆的C，最後求解成功了，那麼表現出來也就是動量不守恆的，如果胡亂定義了一個約束還要求動量守恆其實沒有討論的必要。

約束投影

這塊公式相當多，
首先對於一個等式約束我們要求的就是這個狀態代入約束函數後輸出一個0
即

\[C(p+\Delta p) = 0 \]

而這個函數可以一階泰勒展開得到

\[C(p+\Delta p) \approx C(p) + \nabla_p C(p) \cdot \Delta p \]

我們要把delta P限制在約束函數的梯度方向即

\[\Delta p = \lambda \nabla_p C(p) \]

這個式子代入上一個式子得到λ
得到

\[\Delta p = - \frac{C(p)}{|\nabla_pC(p)|^2} \nabla_p C(p) \]

這就是對於單個約束函數，整個系統接下來要做的改變。
而對於單個點p_i來說，前面的係數λ是一樣的，後面的梯度方向則變成了這個delta P在單個點的分量，其實是梯度的某幾個維度，也就是關於pi的偏導數
前面的係數用s表示，

\[\Delta p_i = - s \nabla_{p_i} C(p) \]

\[ s = \frac{C(p)}{\sum_j {|\nabla_{p_j}C(p) |^2}} \]

其實這裡下面的分母就是梯度的模長，梯度的模長等於各個偏導數分量的平方和。
特別注意下維度：
這裡的p是指所有的點的位置，維度是點的數量 n * 3 ，其中delta P會保證和約束函數的梯度一樣的方向，約束函數的梯度也是一個 n * 3的向量。

delta P1則表示delta P向量的第一點分量，是一個3分量的向量，其實就是第一個點的位移，他應該保證和約束函數在這個第一個點的偏導數保證方向相同。

簡單約束舉例

這裡就不舉例了，可以檢視參考文章裡面的距離約束
https://zhuanlan.zhihu.com/p/48737753

剛度體現

公式裡還有一個k引數，這個引數的作用意會一下就是，我們有一條約束方程，要以怎麼樣的速度去滿足這條方程，如果k等於1則說明直接滿足，不留餘地，等於0-1的值，則說明可能是跑幾輪才能比較接近精確解。
剛度越大，越快滿足約束方程，也就和剛度體現的越硬對應了。
一個問題是，他自己指出材料的剛度依賴時間步長，定步長就不太會有這些問題，verlet積分的推斷也是依賴定步長。
多次迭代下誤差為

\[ \Delta p (1-k)^{n_s} \]

另外，為了消除多次迭代中殘差隨著次數指數級增長的問題，他這裡用的一個非常怪的式子反推了，然後誤差變為線性依賴於k而不依賴於ns,ns是迭代次數。

\[k' = 1-(1-k)^{\frac{1}{n_s}} \]

最終誤差來到

\[\Delta p (1-k')^{n_s} = \Delta p(1-k) \]

遊戲模擬——Position based dynamics