第二講：動力學方程的數(shù)值積分

計算機動畫的本質(zhì)上是解微分方程，因此本文主要介紹微分方程常用的一些解法。

particle system

首先從最簡單的內(nèi)容開始，也就是particle system。

上面是particle system的一個簡單的案例，也就是一個彈簧+一個block，block在彈簧的作用下左右運動，其中彈簧不受力的情況下block所在的位置為原點，這里用q表示位置，那么上述狀態(tài)可以用q=0來表示，定義向左為負，向右為正，相應的速度與加速度也可以基于下面的公式得到：

\dot{q} = \frac{dq}{dt}
\ddot{q} = \frac{d\dot{q}}{dt}

動能kinetic energyE_k與勢能potential energyE_p則由如下的公式進行表示：

E_k = \frac{1}{2}m\dot{q}^2
E_p = \frac{1}{2}kq^2

其中m與k分別表示block的質(zhì)量與彈簧的彈性系數(shù)，而彈簧的彈力可通過下述公式計算得到：
f = -kq

根據(jù)前面的公式可以知道，彈力是conservative force的一種。

所謂的conservative force指的是將一個particle從一個點移動到另一個點做的功是跟移動的路徑無關(guān)的力，在這個定義下，如果將一個particle沿著任意路徑移動最后回到原點，那么力的整體做功（微分下可以用\int f(x) * dx表示）為0，在conservative force的作用下，不論沿著哪條路徑，只要保證首尾點是相同的，那么都會導致勢能的變化幅度是相同的，常見的力中，重力跟彈力是conservative force，而摩擦力則不是conservative force，從這個角度來理解，conservative force是跟位置相關(guān)的力，而這種力作用的效果就是將位置的變化轉(zhuǎn)換為勢能的變化。

一個力場（force field）如果滿足下面三個公式中的任意一個，那么就可以說明這個力場是conservative的：

  1. 力的curl是0向量
    \vec{\nabla} \times \vec{F} = 0
    在2D空間中，上述公式可以簡化為：
    \frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y} = 0

  2. 閉環(huán)移動時，做功為0
    W = \oint \vec{F} \mathrmu0z1t8os \vec{r} = 0

  3. 力可以通過對勢能的梯度取反來表示：
    \vec{F} = -\vec{\nabla}\Phi

前面兩個公式計算比較復雜，通常我們會使用第三個公式進行conservative force的判斷，根據(jù)前面彈簧粒子的公式，我們可以判斷彈力就是conservative force。

再說回彈簧粒子系統(tǒng)，根據(jù)能量守恒定律，我們有：
\frac{1}{2}m\dot{q}^2 + \frac{1}{2}kq^2 = C
其中C是常量，對這個公式的左右兩邊對時間進行求導，就可以得到：
m\dot{q}\ddot{q} + kq\dot{q}=\dot{q}(m\ddot{q}+kq) = 0
上述這個公式有兩個解：
\dot{q} = 0 \\ m\ddot{q} = - kq
其中第一個解的含義是速度\dot{q}不變，也就是物件靜止，這種情況我們是不關(guān)心的，我們主要關(guān)心第二個解，這個解我們換個寫法：
\ddot{q} = \frac{- kq}{m}
用中文的話翻譯一下，就是物體運動的加速度大小與物體的受力成正比，跟物體的質(zhì)量成反比，這就是牛頓第二定律（作為回顧：牛頓第一定律也稱為慣性定律，可以表述為任何物體都要保持勻速直線運動或靜止狀態(tài)，直到外力迫使它改變運動狀態(tài)為止；牛頓第三定律可以表述為相互作用的兩個物體之間的作用力和反作用力總是大小相等，方向相反，作用在同一條直線上）。

下面再來看一個復雜一點的案例，block從一個變成兩個：

動能公式就變成了：
E_k = \frac{1}{2}m_1\dot{q_1}^2 + \frac{1}{2}m_2\dot{q_2}^2

用矩陣來表示就可以寫成：
E_k = \frac{1}{2} \left[ \begin{matrix} \dot{q_1} \\ \dot{q_2} \end{matrix} \right]^T \left[ \begin{matrix} m_1 & 0 \\ 0 & m_2 \end{matrix} \right] \left[ \begin{matrix} \dot{q_1} \\ \dot{q_2} \end{matrix} \right]

這里需要注意的是，我們前面說過，單個block的情況下，我們是以彈簧原始長度時block所在的位置為原點的，而實際上我們這里選擇的坐標系是比較隨性的，而不是有固定的一套標準，而使用不同的坐標系，得到的公式是不一樣的，雖然他們輸出的結(jié)果是一樣的（能量不隨坐標系的變化），而這給了我們一個啟發(fā)，在實際使用中要謹慎選擇坐標系，不同的坐標系計算公式不一樣，導致計算效率會受影響，而這對于結(jié)果的輸出并沒有任何收益。

數(shù)值積分

重寫一下前面的牛頓第二定律的公式：

\ddot{q} = \frac{- kq}{m}

這是一個以q為變量的二階微分方程，在物理動畫模擬中，我們最終需要的是計算各個時刻的位置，也就是解出每個time step下的q，而求解微分方程的方式就是數(shù)值積分。

為了解前面這個二階微分方程，我們這里嘗試將之改寫成一個一階的常微分方程（目的是降低求解的復雜度？），令
y = \left[ \begin{matrix} q \\ \dot{q} \end{matrix} \right]
那么有
\frac{\textu0z1t8osy}{\textu0z1t8ost} = \left[ \begin{matrix} \dot{q} \\ \ddot{q} \end{matrix} \right] = \left[ \begin{matrix} \dot{q} \\ \frac{-kq}{m} = \frac{F(q, \dot{q})}{m} \end{matrix} \right] = \mathscr{F} (y)
上面為了進行一般化處理，用F(q, \dot{q})替代了-kq，那么怎么來理解這個y以及y的微分呢？

如下圖所示，我們以q作為橫坐標，\dot{q}作為縱坐標繪制一個坐標系（這個坐標系這里稱之為phase space，當然，這里只是一個近似，真實的phase space的縱坐標是動量m\dot{q}，不過不影響），當我們?nèi)↑cy = (0, \dot{q})來看，其微分結(jié)果為(\dot{q}, 0)，也就是一條水平方向（向右）的向量，如下面所示：

同樣，如果我們?nèi)↑cy = (q, 0)來看，其微分結(jié)果為(0, \frac{-kq}{m})，也就是一條垂直方向（向下）的向量，如下面所示：

實際上，對單block的彈簧系統(tǒng)而言，對于每個y坐標，在phase space中得到的y的微分，稱為y的速度向量，都可以得到這個向量是沿著順時針方向的，如果我們做特殊處理，取k = m，那么將所有點連起來，我們就能得到一個正圓，如下圖所示：

不同半徑的正圓對應于彈簧和諧運動的振幅，通過這個正圓也可以看到，在q為0的時候，速度的絕對值最大，反之速度為0的時候，則是位置的絕對值最大的時候；另外，對于任意一個y，其速度\frac{\mathrmu0z1t8osy}{\mathrmu0z1t8ost}都是這個正圓上的切線，也就是說下一刻的位置同樣是在正圓上，如果沒有摩擦力的存在，那么就意味著這個block將在這個圓上永無止盡的運動下去。

在進行數(shù)值積分方法介紹之前，我們先來看下不同的積分方法的評價標準，總的來說，對于一個數(shù)值積分方法好不好，我們通常有如下三個評判標準：

  1. stability，穩(wěn)定性，整個系統(tǒng)會隨著時間的變遷變得穩(wěn)定，不論是勻速運動，還是隨著能量的衰減逐漸趨于靜止，或者在某個有限的能量的范圍內(nèi)來回振蕩
  2. accuracy，精確度，指的是我們通過數(shù)值積分模擬出來的結(jié)果跟真實的精確解之間的誤差
  3. convergence，收斂度，指的是我們在進行數(shù)值積分模擬計算的時候，模擬的time step在不斷收縮的情況下，其模擬的精度是否不斷逼近真實解。

將上面的公式做一下近似，這里用explicit euler（顯式歐拉）公式進行模擬，即\scr{F}(y)中y取k：
\frac{y(k+1) - y(k)}{h} \approx \frac{\textu0z1t8osy}{\textu0z1t8ost} = \mathscr{F} (y(k))
從而得到
y(k+1) = h \cdot \mathscr{F}(y(k)) + y(k)

按照上面這個公式，由于\scr{F}(y)是在y點的切線，也就是說計算得到的下一個坐標點將落在當前坐標點所規(guī)范的圓形外面，那么經(jīng)過這個公式計算出來的運動軌跡將不能保證是圓形的，而是半徑不斷擴大的螺旋曲線，也就是說，這個方法是不穩(wěn)定的，即系統(tǒng)的能量越來越大（半徑對應著振幅）：

在顯式歐拉方法下，我們看下整體方案的誤差，根據(jù)泰勒展式：
y(t_0+h) = y(t_0) + \dot y(t_0) \cdot h + \ddot y(t_0) \cdot \frac{h^2}{2} + ...
可以看到，顯式歐拉只使用了\ddot y之前的部分，也就是說是屬于一階微分方程，其誤差自然就可以算出來等于O(h^2)了。

如果前面不用顯式歐拉，而是使用implicit euler（隱式歐拉），即\scr{F}(y)中y取k+1：
\frac{y(k+1) - y(k)}{h} \approx \frac{\textu0z1t8osy}{\textu0z1t8ost} = \mathscr {F}(y(k+1))

在這種情況下，右邊是一個方程，在彈簧系統(tǒng)中，這個方程是可以顯式表達的，但是在其他的系統(tǒng)尤其是復雜系統(tǒng)中，這個方程是沒有辦法表達出來的，也就是說這個公式通過解析方法是沒有辦法計算得到的。

y(k+1) - h\cdot \mathscr{F}(y(k+1)) = y(k)

根據(jù)上面的公式，仿造前面顯式歐拉的邏輯，我們可以看出來，y(k)是比y(k+1)半徑要大的，也就是說，在這個系統(tǒng)中，半徑在不斷減小，能量在不斷衰退，最終衰減到0，那么就處于靜止，也就是說這種方式得到的是一個穩(wěn)定的解（當然，這里需要說清楚的是，不是所有的系統(tǒng)使用隱式歐拉都是穩(wěn)定的，只有那些能量守恒的系統(tǒng)才是穩(wěn)定的），不過呢由于\scr{F}沒有解析解，因此求解會很麻煩。

下面要介紹的是一個叫做symplectic euler（辛歐拉，或者半顯式歐拉，半隱式歐拉方案）的數(shù)值積分方法，根據(jù)前面的公式，我們有：
y = \left[ \begin{matrix} q \\ \dot{q} \end{matrix} \right]
\dot{y} = \left[ \begin{matrix} \dot{q} \\ \ddot{q} \end{matrix} \right] = \scr{F}(q, \dot{q})
將之轉(zhuǎn)換為微分方程：
\frac{1}{h} \left[ \begin{matrix} q(k+1) - q(k) \\ \dot{q}(k+1) - \dot{q}(k) \end{matrix} \right] = \left[ \begin{matrix} \dot{q} \\ \frac{F}{m}(q, \dot{q}) \end{matrix} \right]
又回到前面的問題，右邊的向量中的q是應該取k還是k+1時刻的值，前面我們說過，取k的值是顯式歐拉的方案，好處是數(shù)值是已知的，這個方程可以解出來的，缺點則是系統(tǒng)不穩(wěn)定，而取k+1的時候則是隱式歐拉的方案，好處是方案是穩(wěn)定的，而缺點則是方程求解困難，基于這一點，我們這里考慮使用一個混合的方案，即上述方程中第二個分量使用顯式歐拉方案，第一個分量則使用隱式歐拉方案，如下所示：
\frac{1}{h} \left[ \begin{matrix} q(k+1) - q(k) \\ \dot{q}(k+1) - \dot{q}(k) \end{matrix} \right] = \left[ \begin{matrix} \dot{q}(k+1) \\ \frac{F}{m}(q(k), \dot{q}(k)) \end{matrix} \right]
根據(jù)上面的公式，我們可以先根據(jù)第二個分量的公式，計算出\dot{q}(k+1)，之后這個數(shù)值可以代入到第一個公式，求得q(k+1)，這里借用了顯式歐拉的好處，只是需要確定的是，是否同時帶入了顯式歐拉的缺點——不穩(wěn)定。

在彈簧系統(tǒng)中，上述公式可以表示為：
\frac{1}{h} \left[ \begin{matrix} q(k+1) - q(k) \\ \dot{q}(k+1) - \dot{q}(k) \end{matrix} \right] = \left[ \begin{matrix} \dot{q}(k+1) \\ \frac{-k\cdot q(k)}{m} \end{matrix} \right]
從而求得：
\dot{q}(k+1) = \dot{q}(k) - \frac{h \cdot k \cdot q(k)}{m} \\ q(k+1) = q(k) + h \cdot \dot{q}(k+1) = q(k) + \dot{q}(k) \cdot h- \frac{ k \cdot q(k)}{m} h^2
根據(jù)這公式通過不斷迭代可以得到：
q(k+1) = q(k) - \frac{kh^2}{m}\sum_{i=0}^k q(i) +\dot{q}(0) \cdot h
好像也不能說明啥，寫成矩陣形式：
\left[ \begin{matrix} q(k+1) \\ \dot{q}(k+1) \end{matrix} \right] = \left[ \begin{matrix} 1-\frac{k}{m} h^2 && h \\ -\frac{kh}{m} && 1 \end{matrix} \right] \left[ \begin{matrix} q(k) \\ \dot{q}(k) \end{matrix} \right]
這個等式右邊的2x2的矩陣，我們可以很容易計算出其行列式=1，而根據(jù)線性代數(shù)的理論，我們可以知道，矩陣的行列式恰好等于變換前后兩個形狀的面積之比，也就是說，某個區(qū)域的數(shù)據(jù)集合經(jīng)過這個矩陣變化后得到的新的區(qū)域，兩個區(qū)域的面積之比就恰好等于矩陣的行列式，而這里行列式為1就說明經(jīng)過這個矩陣的變換，數(shù)據(jù)是基本維持穩(wěn)定的，不會膨脹，也不會坍縮，而是能量守恒的（自然也是穩(wěn)定的系統(tǒng)）。

而雖然上述結(jié)論是在彈簧系統(tǒng)中推導出來的，但實際上這個結(jié)論的應用面是可以推廣到一般的（當然，不是所有的系統(tǒng)都可以使用），對于系統(tǒng)內(nèi)僅包含動能和勢能的系統(tǒng)而言，在經(jīng)典力學范疇中，這個結(jié)論是成立的，即通過辛歐拉算法得到的解都是穩(wěn)定的（甚至能量守恒的）。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容